Load libraries and functions

 # install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE) 
start_time <- Sys.time()

library(INLA)          # For fitting integrated nested Laplace approximation (INLA) models
library(tidyverse)     # For data manipulation and visualization
library(sp)            # For spatial data handling
library(spdep)         # For spatial dependence analysis
library(flextable)     # For creating flexible tables
library(RColorBrewer)  # For color palettes
library(rcartocolor)   # For additional color palettes
library(cowplot)       # For combining plots
library(janitor)       # For cleaning data
library(broom)

# Functions ---------------------------------------------------------------

# Standardize a numeric vector
my_std <- function(x) { (x - mean(x,na.rm=TRUE)) / sd(x,na.rm = TRUE)}

# Fit a BYM2 model using INLA
fit_bym2 <- function(table,
                     outcome_var,
                     interest_var,
                     distribution = "poisson",
                     formula_terms = NULL,
                     prec_param =  c(1, 0.001),
                     phi_param = c(0.5, 0.5)) {

# Define priors for precision and phi
    prior <- list(
  prec = list(
    prior = "pc.prec",
    param = prec_param),
  phi = list(
    prior = "pc",
    param = phi_param)
  )
    
  # Construct the INLA formula based on input parameters
  if (is.null(formula_terms)) {
    inla_formula <- as.formula(paste(outcome_var, "~", interest_var, "+ f(ui, model = 'bym2', graph = g, hyper = prior, adjust.for.con.comp = TRUE,
constr = TRUE, scale.model = TRUE)"))
  } else {
    inla_formula <- as.formula(paste(outcome_var, "~", interest_var, "+", formula_terms, "+ f(ui, model = 'bym2', graph = g, hyper = prior, adjust.for.con.comp = TRUE,
constr = TRUE, scale.model = TRUE)"))
  }
    
    
    
#adjust.for.con.comp = TRUE, scale.model = TRUE:  assigns one intercept to each region in addition to using a sum-to-zero constraint for each connected region. By adding an intercept for each region, we infer that the baseline prevalence is different in the disconnected regions.
#The spatial random effects can be interpreted as the area-specific deviation from the region-specific risk in this case. The intercepts for the disconnected regions need to be explicitly specified in INLA:

  
  # Fit the INLA model with BYM2 prior
fit <- inla(inla_formula,
              data = table,
              family = distribution,
              control.compute = list(dic = TRUE, waic = TRUE,cpo=TRUE,config=TRUE),
              control.predictor = list(compute = TRUE),
              verbose = FALSE)
  
  # Return the INLA model object
  return(fit)
}

# Tidy up INLA model fixed effects results
tidy_inla <- function(fixed_effects, exponentiate = FALSE, sigfig =NULL) {
  # Create a tidy data frame
  tidy_results <- fixed_effects %>%
    dplyr::rename(estimate = 'mean',
                  std.error = 'sd',
                  lower_ci = '0.025quant',
                  upper_ci = '0.975quant')
  
  if (exponentiate) {
    tidy_results <- tidy_results %>%
      mutate(estimate = exp(estimate),
             lower_ci = exp(lower_ci),
             upper_ci = exp(upper_ci))
  }
  
  if (!is.null(sigfig)) {
    tidy_results <- tidy_results %>%
      mutate(estimate = signif(estimate, sigfig),
             std.error = signif(std.error, sigfig),
             lower_ci = signif(lower_ci, sigfig),
             upper_ci = signif(upper_ci, sigfig))
  }
  
  return(tidy_results)
}

# prs <- function(x,n.dig = 2){
#   
#   formatC(round(x,n.dig), format='f', big.mark = ",",digits = n.dig)
#   
# }

cma_models <- function(cma_data, cma_name, outcome, deprivation_indicator,distribution = "poisson") {
  
  ###### Spatial Models #####  
  
  ### 1: Unadjusted
  
  print("Fit unadjusted")
  
  m1 <- fit_bym2(cma_data, outcome, deprivation_indicator,distribution)
  m1_fe <- tidy_inla(m1$summary.fixed, exponentiate = TRUE, sigfig =3) %>%
    mutate(
      `IRR (95% CI)` = paste0(estimate, " (", lower_ci, ", ", upper_ci, ")"),
      variable = row.names(.)
    ) %>%
    select(variable, `IRR (95% CI)`)
  
    print("Fit adjusted by road length")

  ### 2: Adjusted by Road Length
  m2 <- fit_bym2(cma_data, outcome, deprivation_indicator, formula_terms = "ln_roads_km_c",distribution)
  m2_fe <- tidy_inla(m2$summary.fixed, exponentiate = TRUE, sigfig =3) %>%
    mutate(
      `IRR (95% CI)` = paste0(estimate, " (", lower_ci, ", ", upper_ci, ")"),
      variable = row.names(.)
    ) %>%
    select(variable, `IRR (95% CI)`)
  
  ### 3: Adjusted by Road Length and other variables
  m3 <- fit_bym2(cma_data, outcome, deprivation_indicator, 
                  formula_terms = "ln_roads_km_c + roads_prop_highway_arterial_z + ale_index_z + canbics_index_z + population_100_c",distribution)
  m3_fe <- tidy_inla(m3$summary.fixed, exponentiate = TRUE, sigfig =3) %>%
    mutate(
      `IRR (95% CI)` = paste0(estimate, " (", lower_ci, ", ", upper_ci, ")"),
      variable = row.names(.)
    ) %>%
    select(variable, `IRR (95% CI)`)

  l <- list(m1, m2, m3)
  names(l) <- c("Unadjusted","Adjusted1","Adjusted2")
  
      print("Fit adjusted by road length and covariates")


  # Compile results into a tibble
  tibble(
    cma = cma_name,
    model_type = c("unadjusted", "adjusted_minimal", "adjusted_full"),
    models = l,
    fixed_effects = list(m1_fe, m2_fe, m3_fe)
  )
         
}


# Clean and format fixed effects results from models
clean_fixed_effects <- function(cma_model_object){
  
  cma_model_object %>% 
    bind_rows(,.id="model_type") %>% 
    mutate(model_type = case_when(model_type == "1"~"Unadjusted IRR (95% CI)",
                                  model_type == "2"~"Minimally Adjusted IRR (95% CI)",
                                  model_type == "3"~"Fully Adjusted IRR (95% CI)",
    )) %>% 
    pivot_wider(id_cols = variable,names_from = model_type,values_from = `IRR (95% CI)` ) %>% 
    mutate(variable = case_when(variable == "ln_roads_km_c" ~"Log(Total KMs of Road)",
                                variable == "roads_prop_highway_arterial_z" ~ "% of roads classifed as Arterial/Highway (SD)",
                                variable == "ale_index_z" ~ "Canadian Active Living Enrivonment Index (SD)",
                                variable == "canbics_index_z" ~ "Canadian Bikeway Safety and Comfort Index (SD)",
                                variable == "population_100_c" ~ "Population (100)",
                                TRUE ~ variable
                                
    )) 
  
}


assign_nearest_neighbors <- function(nb_object, sp_object,make_symmetrical=TRUE) {
    # Calculate centroids for isolated polygons
  centroids <- coordinates(sp_object)  # If using sp
  
  centroids_nn <- knearneigh(centroids, k = 1)

  # Get the list of neighbors
  zero_nb_logical <- sapply(nb_object, function(x) x[1] ==0)
  
  zero_nb_indices <- which(zero_nb_logical)
  
  zero_nn <- centroids_nn$nn[zero_nb_indices, ]
  
  for(i in 1:length(zero_nn)){
    
    nb_object[zero_nb_indices[i]] <- zero_nn[i]
  }
  
  if(make_symmetrical==TRUE){
    
    nb_object <- make.sym.nb(nb_object)
  }
  
  
  return(nb_object)
}

create_descriptive_table <- function(sf_obj, selected_vars) {
  sf_obj %>%
    # Drop geometry and select only the specified variables
    st_drop_geometry() %>%
    select(all_of(selected_vars)) %>%
    # Gather descriptive statistics for each variable
    summarise(across(everything(), list(
      sum = ~sum(.x, na.rm = TRUE),
      mean = ~mean(.x, na.rm = TRUE),
      sd = ~sd(.x, na.rm = TRUE),
      min = ~min(.x, na.rm = TRUE),
      max = ~max(.x, na.rm = TRUE),
      n = ~sum(!is.na(.x)|is.na(.x)),       # Count non-NA values
      missing = ~sum(is.na(.x)),
      n_zero = ~sum(.x==0,na.rm=TRUE)# Count 0 values
    ), .names = "{.col}-{.fn}")) %>%
    # Pivot longer to have each variable as its own row
    pivot_longer(everything(), 
                 names_to = c("variable", ".value"),
                 names_sep = "-") %>%
    # Add a column to identify the sf object
    # Reorder columns for better readability
    select(variable, n, sum, mean, sd, min, max, , missing,n_zero)
}

Load Data and Clean

# Load and clean data -----------------------------------------------------

cma_name <- cancensus::get_census("CA21",
                                  level = "CMA",
                                  regions = list(PR=59)) %>% 
  rename(CMA_UID = GeoUID,
         `CMA Name` = `Region Name`,
         cma_population = Population) %>% 
  arrange(desc(cma_population)) %>% 
  select(CMA_UID, Type, `CMA Name`,cma_population) %>% 
  mutate(`CMA Name` = str_remove_all(`CMA Name`," \\(B\\)|\\(K\\)|\\(D\\)")) %>%
  clean_names() %>% 
  mutate(cma_name = stringr::str_trim(cma_name,"right"))
## Downloading: 920 B     Downloading: 920 B     Downloading: 920 B     Downloading: 920 B     Downloading: 920 B     Downloading: 920 B
# Read spatial data for dissemination areas in British Columbia
bc_da <- st_read(
  "C:/Users/micha/Documents/GitHub/Area-Level-Deprivation-Traffic-Injury/Processed Data/da_v3_2021.gpkg"
) %>%
  # Join the census data with the spatial data
  left_join(cma_name, by = join_by(cma_uid)) %>%
  filter(!is.na(cma_name)) %>% 
  mutate(
    population_100 = population/100,
    
    population_100_c = scale(population_100, scale = FALSE),    # Convert road lengths from meters to kilometers
    total_roads_km = total_roads_m / 1000,
    
    n_casualty_claims_rate = n_casualty_claims/total_roads_km*10,
    n_pedestrian_casualty_claims_rate = n_pedestrian_casualty_claims/total_roads_km*10,
    n_cyclist_casualty_claims_rate = n_cyclist_casualty_claims/total_roads_km*10,

    # Calculate the proportion of roads that are highways or arterials
    roads_prop_highway_arterial = (highway_m + arterial_m) / total_roads_m,
    # Standardize the proportion of highways and arterials
    roads_prop_highway_arterial_z = my_std(roads_prop_highway_arterial),
    # Log-transform the total length of roads in kilometers
    ln_roads_km = ifelse(total_roads_km ==0, NA,log(total_roads_km)),
    ln_roads_km_c = scale(ln_roads_km, scale = FALSE),
    # Calculate the casualty claims rate by total road kilometers
    casualty_claims_rate = n_casualty_claims / total_roads_km * 10,
    # Calculate the cyclist casualty claims rate by total road kilometers
    cyclist_casualty_claims_rate = n_cyclist_casualty_claims / total_roads_km * 10,
    # Calculate the pedestrian casualty claims rate by total road kilometers
    pedestrian_casualty_claims_rate = n_pedestrian_casualty_claims / total_roads_km * 10,
    # Convert journey-to-work total from meters to kilometers
    jtw_total_100 = jtw_total / 100,
    # Standardize the journey-to-work total
    jtw_total_z = my_std(jtw_total),
    # Calculate the proportion of journey-to-work trips that are motor vehicle-related
    jtw_prop_mv = ((jtw_driver + jtw_passenger) / jtw_total),
    # Calculate the proportion of non-motor vehicle journey-to-work trips
    jtw_prop_non_mv = (1 - jtw_prop_mv),
    # Standardize the proportion of non-motor vehicle trips
    jtw_prop_non_mv_z = my_std(jtw_prop_non_mv),
    # Adjust and scale the remote index for visualization
    remote_index_2 = (1 - remote_index) * 100,
    # Log-transform the adjusted remote index
    ln_remote_index = log(remote_index_2),
    # Standardize the VANDIX index
    vandix_z = my_std(vandix),
    # Convert VANDIX quintile to descriptive categories
    vandix_quintile = vandix_quintile %>% as.character(),
    vandix_quintile = case_when(
      vandix_quintile == "1" ~ "1 - Least Deprived",
      vandix_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ vandix_quintile
    ),
    # Convert educational attainment quintile to descriptive categories
    no_highschool_quintile = no_highschool_quintile %>% as.character(),
    no_highschool_quintile = case_when(
      no_highschool_quintile == "1" ~ "1 - Least Deprived",
      no_highschool_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ no_highschool_quintile
    ),
    # Convert unemployment rate quintile to descriptive categories
    unemployment_rate_quintile = unemployment_rate_quintile %>% as.character(),
    unemployment_rate_quintile = case_when(
      unemployment_rate_quintile == "1" ~ "1 - Least Deprived",
      unemployment_rate_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ unemployment_rate_quintile
    ),
    # Convert university degree quintile to descriptive categories
    university_degree_quintile = university_degree_quintile %>% as.character(),
    university_degree_quintile = case_when(
      university_degree_quintile == "1" ~ "1 - Least Deprived",
      university_degree_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ university_degree_quintile
    ),
    # Convert lone-parent family quintile to descriptive categories
    lone_parent_fam_quintile = lone_parent_fam_quintile %>% as.character(),
    lone_parent_fam_quintile = case_when(
      lone_parent_fam_quintile == "1" ~ "1 - Least Deprived",
      lone_parent_fam_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ lone_parent_fam_quintile
    ),
    # Convert average household income quintile to descriptive categories
    hh_avg_income_quintile = hh_avg_income_quintile %>% as.character(),
    hh_avg_income_quintile = case_when(
      hh_avg_income_quintile == "1" ~ "1 - Least Deprived",
      hh_avg_income_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ hh_avg_income_quintile
    ),
    # Convert home ownership quintile to descriptive categories
    home_owner_quintile = home_owner_quintile %>% as.character(),
    home_owner_quintile = case_when(
      home_owner_quintile == "1" ~ "1 - Least Deprived",
      home_owner_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ home_owner_quintile
    ),
    # Convert participation rate quintile to descriptive categories
    participation_rate_quintile = participation_rate_quintile %>% as.character(),
    participation_rate_quintile = case_when(
      participation_rate_quintile == "1" ~ "1 - Least Deprived",
      participation_rate_quintile == "5" ~ "5 - Most Deprived",
      TRUE ~ participation_rate_quintile
    ),
    # Convert CANBICS class to descriptive categories
    canbics_class_c = case_when(
      canbics_class == 1 | canbics_class == "2" ~ "1 - Low",
      canbics_class == 5 | canbics_class == "4" ~ "3 - High",
      canbics_class == 3 ~ "2 - Medium"
    ),
    # Standardize the CANBICS index
    canbics_index_z = my_std(canbics_index),
    
    # Convert ALE class to descriptive categories
    ale_class_c = case_when(
      ale_class == 1 | ale_class == "2" ~ "1 - Low",
      ale_class == 5 | ale_class == "4" ~ "3 - High",
      ale_class == 3 ~ "2 - Medium"
    ),
    # Standardize the ALE index
    ale_index_z = my_std(ale_index),
    # Calculate VanDIX 
    z_no_highschool_prevalance_w = scale(no_highschool_prevalance),
    z_university_degree_prevalance_w = scale(university_degree_prevalance * -1),
    z_unemployment_rate_w = scale(unemployment_rate),
    z_lone_parent_fam_prevalence_w = scale(lone_parent_fam_prevalence),
    z_hh_avg_income_w = scale(hh_avg_income * -1),
    z_home_owner_prevalence_w = scale(home_owner_prevalence * -1),
    z_participation_rate_w = scale(participation_rate * -1),
    # Compute VanDIX (Area Level Deprivation Index) score
    vandix = as.numeric(
      z_hh_avg_income_w * 0.089 +
        z_home_owner_prevalence_w * 0.089 +
        z_lone_parent_fam_prevalence_w * 0.143 +
        z_no_highschool_prevalance_w * 0.25 +
        z_university_degree_prevalance_w * 0.179 +
        z_unemployment_rate_w * 0.214 +
        z_participation_rate_w * 0.036
    ),
    vandix_z_c = case_when(vandix_z< -5 ~ "<-5",
                                    vandix_z>= -5 & vandix_z< -4 ~"-5 - -4",
                                    vandix_z>= -4 & vandix_z < -3 ~"-4 - -3",
                                    vandix_z>= -3 & vandix_z < -2 ~"-3 - -2",
                                    vandix_z>= -2 & vandix_z < -1 ~"-2 - -1",
                                    vandix_z>= -1 & vandix_z < 0 ~"-1 - 0",
                                    vandix_z>= 0 & vandix_z < 1 ~"0 - 1",
                                    vandix_z>= 1 & vandix_z<2 ~"1 - 2",
                                    vandix_z>= 2 & vandix_z<3 ~"2 - 3",
                                    vandix_z>= 3 & vandix_z<4 ~"3 - 4",
                                    vandix_z>= 4 & vandix_z<5 ~"4 - 5",
                                    vandix_z >= 5 ~ ">5",
      ),
      vandix_z_c  = factor(vandix_z_c,levels = c("<-5","-5 - -4","-4 - -3" , "-3 - -2" , "-2 - -1" , "-1 - 0",  "0 - 1" ,  "1 - 2",  "2 - 3", "3 - 4","4 - 5",">5"))
  ) %>% 
    mutate(region_name = case_when(
    cma_name %in% c("Prince Rupert", "Terrace") ~ "Northwest",
    cma_name %in% c("Fort St. John", "Dawson Creek", "Prince George", "Quesnel", "Williams Lake") ~ "North Central",
    cma_name %in% c("Kamloops", "Salmon Arm") ~ "Kamloops-Salmon Arm",
    cma_name %in% c("Vernon", "Kelowna", "Penticton") ~ "Okanagan",
    cma_name %in% c("Cranbrook", "Nelson", "Trail") ~ "Southeast",
    cma_name %in% c("Vancouver", "Squamish") ~ "Vancouver-Squamish",
    cma_name %in% c("Abbotsford - Mission", "Chilliwack") ~ "Fraser Valley",
    cma_name %in% c("Campbell River", "Courtenay", "Port Alberni", "Parksville", "Nanaimo", "Ladysmith", "Duncan","Powell River") ~ "Central Island-Powell River",
    cma_name == "Victoria" ~ "Victoria",
    TRUE ~ NA_character_  # Handle cases where cma_name doesn't match any of the specified values
  ))
## Reading layer `da_v3_2021' from data source 
##   `C:\Users\micha\Documents\GitHub\Area-Level-Deprivation-Traffic-Injury\Processed Data\da_v3_2021.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 7848 features and 130 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 273876.1 ymin: 367780.7 xmax: 1870608 ymax: 1735671
## Projected CRS: NAD83 / BC Albers
#split data into list of CMAs
bc_cma <- bc_da %>% 
  split(.$cma_name)

bc_region <- bc_da %>% 
  split(.$region_name)

cma_pop <- bc_da %>%
  st_drop_geometry() %>% 
  group_by(cma_name) %>% 
  summarise(n_da = n(),
            population = sum(population,na.rm=TRUE)) %>% 
  arrange(desc(population))

region_pop <- bc_da %>%
  st_drop_geometry() %>% 
  group_by(region_name) %>% 
  summarise(n_da = n(),
            population = sum(population,na.rm=TRUE)) %>% 
  arrange(desc(population))



region_descriptives <- map(
  bc_region,
  create_descriptive_table,
  c(
    "n_claims",
    "n_casualty_claims",
    "n_cyclist_claims",
    "n_cyclist_casualty_claims",
    "n_pedestrian_claims",
    "n_pedestrian_casualty_claims",
    "population",
    "total_roads_km",
    "roads_prop_highway_arterial",
    "no_highschool_prevalance",
    "unemployment_rate",
    "hh_avg_income",
    "participation_rate",
    "university_degree_prevalance",
    "lone_parent_fam_prevalence",
    "home_owner_prevalence",
    "vandix",
    "roads_prop_highway_arterial",
    "ale_index",
    "canbics_index"
  )
) %>% 
  bind_rows(.id = "region") %>% 
  arrange(desc(n)) %>% 
  mutate(p_zero = n_zero/n*100)


cma_descriptives <- map(
  bc_cma,
  create_descriptive_table,
  c(
    "n_claims",
    "n_casualty_claims",
    "n_cyclist_claims",
    "n_cyclist_casualty_claims",
    "n_pedestrian_claims",
    "n_pedestrian_casualty_claims",
    "population",
    "total_roads_km",
    "roads_prop_highway_arterial",
    "no_highschool_prevalance",
    "unemployment_rate",
    "hh_avg_income",
    "participation_rate",
    "university_degree_prevalance",
    "lone_parent_fam_prevalence",
    "home_owner_prevalence",
    "vandix",
    "roads_prop_highway_arterial",
    "ale_index",
    "canbics_index"
  )
) %>% 
  bind_rows(.id = "cma") %>% 
  arrange(desc(n)) %>% 
  mutate(p_zero = n_zero/n*100)

Spatial Models

Vancouver - Squamish

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Vancouver-Squamish

n_claims

3,622

693,530.0

191.5

363.2

0.0

6,315.0

0

12

0.3

n_casualty_claims

149,692.0

41.3

84.8

0.0

1,896.0

0

133

3.7

n_cyclist_claims

7,656.0

2.1

4.7

0.0

115.0

0

1,457

40.2

n_cyclist_casualty_claims

4,968.0

1.4

3.2

0.0

72.0

0

1,859

51.3

n_pedestrian_claims

9,539.0

2.6

5.1

0.0

80.0

0

1,290

35.6

n_pedestrian_casualty_claims

7,502.0

2.1

4.0

0.0

68.0

0

1,458

40.3

population

2,667,057.0

736.6

532.3

0.0

8,800.0

1

10

0.3

total_roads_km

15,056.8

4.2

5.7

0.0

176.7

0

2

0.1

roads_prop_highway_arterial

577.0

0.2

0.2

0.0

1.0

2

1,426

39.4

no_highschool_prevalance

476.5

0.1

0.1

0.0

0.8

300

16

0.4

unemployment_rate

19,517.2

5.9

3.6

0.0

40.0

300

339

9.4

hh_avg_income

279,328,113.0

85,421.4

34,171.3

0.0

560,267.0

352

3

0.1

participation_rate

215,750.1

64.9

10.1

0.0

94.6

300

1

0.0

university_degree_prevalance

1,890.5

0.6

0.1

0.1

1.0

300

0

0.0

lone_parent_fam_prevalence

513.7

0.2

0.1

0.0

0.7

299

15

0.4

home_owner_prevalence

2,229.1

0.7

0.2

0.0

1.3

300

36

1.0

vandix

-473.2

-0.1

0.6

-2.1

3.9

352

0

0.0

ale_index

5,407.4

1.6

3.6

-2.1

25.1

310

13

0.4

canbics_index

15,159.2

4.6

4.1

0.0

24.3

305

72

2.0

Maps

claims <- ggplot() + 
  geom_sf(data = vancouver_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = vancouver_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = vancouver_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

vancouver_sp <- as(vancouver_sf, "Spatial")
vancouver_sp$ui <- 1:nrow(vancouver_sp@data)

coords <- coordinates(vancouver_sp)
vancouver_nb <- poly2nb(vancouver_sp, queen = TRUE)
## Warning in poly2nb(vancouver_sp, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in poly2nb(vancouver_sp, queen = TRUE): neighbour object has 6 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
vancouver_nb
## Neighbour list object:
## Number of regions: 3622 
## Number of nonzero links: 22972 
## Percentage nonzero weights: 0.1751064 
## Average number of links: 6.342352 
## 1 region with no links:
## 2938
## 6 disjoint connected subgraphs
#assign nearest neighbour for no links

vancouver_nb <- assign_nearest_neighbors(vancouver_nb,vancouver_sp)

plot(vancouver_sp, border = grey(0.5))
plot(vancouver_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(vancouver_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(vancouver_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  vancouver_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 24.441, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.318516e-01     -2.761668e-04      9.020552e-05
cyc_cc_mi <- moran.test(vancouver_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  vancouver_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 38.928, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      3.683422e-01     -2.761668e-04      8.966834e-05
pd_cc_mi <- moran.test(vancouver_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  vancouver_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 32.695, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      3.122914e-01     -2.761668e-04      9.139466e-05

BYM2 Models

nb2INLA("vancouver.adj", vancouver_nb)
g <- inla.read.graph(filename = "vancouver.adj")


# Model Set 1: Total Casualty Claim Crashes

van_models_1 <- cma_models(vancouver_sp@data, "Vancouver-Squamish", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
van_all <- clean_fixed_effects(van_models_1$fixed_effects)

map(van_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 18.7, Running = 5.27, Post = 0.528, Total = 24.5 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 2.802 0.016      2.771    2.802      2.834 2.802   0
## vandix_z    0.211 0.033      0.147    0.211      0.276 0.211   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.311 0.014      0.284    0.310      0.339 0.310
## Phi for ui       0.794 0.021      0.751    0.795      0.833 0.797
## 
## Deviance Information Criterion (DIC) ...............: 23817.67
## Deviance Information Criterion (DIC, saturated) ....: 7368.64
## Effective number of parameters .....................: 3420.00
## 
## Watanabe-Akaike information criterion (WAIC) ...: 23712.55
## Effective number of parameters .................: 2288.45
## 
## Marginal log-Likelihood:  -15034.43 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 18.6, Running = 4.79, Post = 0.636, Total = 24 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   3.254 0.015      3.223    3.254      3.284 3.254   0
## vandix_z      0.327 0.026      0.276    0.327      0.378 0.327   0
## ln_roads_km_c 1.294 0.029      1.238    1.294      1.351 1.294   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.508 0.020      0.470    0.508      0.548 0.507
## Phi for ui       0.811 0.016      0.778    0.811      0.841 0.812
## 
## Deviance Information Criterion (DIC) ...............: 23705.93
## Deviance Information Criterion (DIC, saturated) ....: 7256.88
## Effective number of parameters .....................: 3264.92
## 
## Watanabe-Akaike information criterion (WAIC) ...: 23514.43
## Effective number of parameters .................: 2154.47
## 
## Marginal log-Likelihood:  -14241.75 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 18.2, Running = 5.04, Post = 0.736, Total = 24 
## Fixed effects:
##                                mean    sd 0.025quant 0.5quant 0.975quant  mode
## (Intercept)                   2.989 0.017      2.957    2.989      3.022 2.989
## vandix_z                      0.222 0.022      0.180    0.222      0.265 0.222
## ln_roads_km_c                 0.990 0.027      0.937    0.990      1.043 0.990
## roads_prop_highway_arterial_z 0.593 0.016      0.562    0.593      0.624 0.593
## ale_index_z                   0.109 0.037      0.036    0.109      0.182 0.109
## canbics_index_z               0.209 0.037      0.137    0.209      0.282 0.209
## population_100_c              0.017 0.003      0.012    0.017      0.022 0.017
##                               kld
## (Intercept)                     0
## vandix_z                        0
## ln_roads_km_c                   0
## roads_prop_highway_arterial_z   0
## ale_index_z                     0
## canbics_index_z                 0
## population_100_c                0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.838 0.038      0.765    0.837      0.915 0.835
## Phi for ui       0.777 0.022      0.733    0.778      0.818 0.780
## 
## Deviance Information Criterion (DIC) ...............: 23560.95
## Deviance Information Criterion (DIC, saturated) ....: 7111.89
## Effective number of parameters .....................: 3108.88
## 
## Watanabe-Akaike information criterion (WAIC) ...: 23296.41
## Effective number of parameters .................: 2021.13
## 
## Marginal log-Likelihood:  -13603.57 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

van_models_2 <- cma_models(vancouver_sp@data,"Vancouver-Squamish","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
van_cyclist <- clean_fixed_effects(van_models_2$fixed_effects)

map(van_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 18.3, Running = 3.99, Post = 0.406, Total = 22.6 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.523 0.032     -0.588   -0.523     -0.461 -0.523   0
## vandix_z     0.148 0.037      0.075    0.148      0.221  0.148   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.567 0.042      0.489    0.565      0.653 0.562
## Phi for ui       0.633 0.048      0.534    0.635      0.724 0.638
## 
## Deviance Information Criterion (DIC) ...............: 9767.72
## Deviance Information Criterion (DIC, saturated) ....: 5213.94
## Effective number of parameters .....................: 1755.44
## 
## Watanabe-Akaike information criterion (WAIC) ...: 9992.02
## Effective number of parameters .................: 1391.74
## 
## Marginal log-Likelihood:  -4164.36 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 18.2, Running = 4.69, Post = 0.463, Total = 23.4 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -0.134 0.026     -0.186   -0.134     -0.083 -0.134   0
## vandix_z       0.217 0.032      0.154    0.217      0.280  0.217   0
## ln_roads_km_c  1.072 0.035      1.003    1.072      1.141  1.072   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.622 0.047      0.532    0.621      0.718 0.621
## Phi for ui       0.934 0.026      0.874    0.938      0.974 0.945
## 
## Deviance Information Criterion (DIC) ...............: 9035.13
## Deviance Information Criterion (DIC, saturated) ....: 4481.34
## Effective number of parameters .....................: 1133.23
## 
## Watanabe-Akaike information criterion (WAIC) ...: 9017.22
## Effective number of parameters .................: 857.84
## 
## Marginal log-Likelihood:  -3785.23 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 18.4, Running = 5.09, Post = 0.709, Total = 24.2 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -0.291 0.031     -0.352   -0.291     -0.231
## vandix_z                       0.185 0.032      0.123    0.185      0.248
## ln_roads_km_c                  0.844 0.039      0.767    0.844      0.921
## roads_prop_highway_arterial_z  0.301 0.024      0.254    0.301      0.348
## ale_index_z                    0.041 0.045     -0.048    0.041      0.130
## canbics_index_z                0.092 0.048     -0.001    0.092      0.186
## population_100_c               0.022 0.003      0.016    0.022      0.028
##                                 mode kld
## (Intercept)                   -0.291   0
## vandix_z                       0.185   0
## ln_roads_km_c                  0.844   0
## roads_prop_highway_arterial_z  0.301   0
## ale_index_z                    0.041   0
## canbics_index_z                0.092   0
## population_100_c               0.022   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.822 0.069      0.692    0.819      0.964 0.817
## Phi for ui       0.871 0.036      0.790    0.874      0.931 0.880
## 
## Deviance Information Criterion (DIC) ...............: 8930.87
## Deviance Information Criterion (DIC, saturated) ....: 4377.08
## Effective number of parameters .....................: 1046.94
## 
## Watanabe-Akaike information criterion (WAIC) ...: 8897.21
## Effective number of parameters .................: 790.56
## 
## Marginal log-Likelihood:  -3706.06 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

van_models_3 <- cma_models(vancouver_sp@data,"Vancouver-Squamish","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
van_pedestrian <- clean_fixed_effects(van_models_3$fixed_effects)

map(van_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 18, Running = 4.19, Post = 0.429, Total = 22.7 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.035 0.027     -0.089   -0.034      0.018 -0.034   0
## vandix_z     0.262 0.036      0.191    0.262      0.333  0.262   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.496 0.035      0.430    0.494      0.567 0.492
## Phi for ui       0.663 0.044      0.573    0.664      0.746 0.666
## 
## Deviance Information Criterion (DIC) ...............: 11940.16
## Deviance Information Criterion (DIC, saturated) ....: 5994.58
## Effective number of parameters .....................: 2208.49
## 
## Watanabe-Akaike information criterion (WAIC) ...: 12274.22
## Effective number of parameters .................: 1755.74
## 
## Marginal log-Likelihood:  -5474.94 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 18.5, Running = 4.56, Post = 0.393, Total = 23.5 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   0.319 0.022      0.275    0.319      0.363 0.319   0
## vandix_z      0.304 0.032      0.242    0.304      0.367 0.304   0
## ln_roads_km_c 0.997 0.034      0.930    0.997      1.065 0.997   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                  mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.45 0.028      0.396    0.450      0.506 0.450
## Phi for ui       0.95 0.019      0.907    0.953      0.979 0.957
## 
## Deviance Information Criterion (DIC) ...............: 11175.53
## Deviance Information Criterion (DIC, saturated) ....: 5229.93
## Effective number of parameters .....................: 1630.40
## 
## Watanabe-Akaike information criterion (WAIC) ...: 11151.04
## Effective number of parameters .................: 1185.45
## 
## Marginal log-Likelihood:  -5141.60 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 18.5, Running = 5.45, Post = 0.451, Total = 24.4 
## Fixed effects:
##                                mean    sd 0.025quant 0.5quant 0.975quant  mode
## (Intercept)                   0.115 0.028      0.060    0.115      0.168 0.115
## vandix_z                      0.268 0.031      0.208    0.268      0.328 0.268
## ln_roads_km_c                 0.717 0.039      0.641    0.717      0.793 0.717
## roads_prop_highway_arterial_z 0.394 0.022      0.351    0.394      0.438 0.394
## ale_index_z                   0.232 0.046      0.141    0.232      0.322 0.232
## canbics_index_z               0.042 0.047     -0.051    0.042      0.135 0.042
## population_100_c              0.027 0.003      0.021    0.027      0.034 0.027
##                               kld
## (Intercept)                     0
## vandix_z                        0
## ln_roads_km_c                   0
## roads_prop_highway_arterial_z   0
## ale_index_z                     0
## canbics_index_z                 0
## population_100_c                0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.713 0.059      0.602    0.712      0.835 0.709
## Phi for ui       0.850 0.039      0.763    0.853      0.917 0.860
## 
## Deviance Information Criterion (DIC) ...............: 11047.95
## Deviance Information Criterion (DIC, saturated) ....: 5102.36
## Effective number of parameters .....................: 1480.57
## 
## Watanabe-Akaike information criterion (WAIC) ...: 10986.17
## Effective number of parameters .................: 1069.75
## 
## Marginal log-Likelihood:  -4953.61 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
van_results <- bind_rows(van_all %>% mutate(Region = "Vancouver-Squamish",Outcome = "All Injury Claims"),
                         van_cyclist %>% mutate(Region = "Vancouver-Squamish",Outcome = "Cyclist Injury Claims"),
                         van_pedestrian %>% mutate(Region = "Vancouver-Squamish",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

Victoria

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Victoria

n_claims

581

67,182.0

115.6

176.1

0.0

1,665.0

0

3

0.5

n_casualty_claims

12,338.0

21.2

34.2

0.0

305.0

0

26

4.5

n_cyclist_claims

1,471.0

2.5

4.4

0.0

53.0

0

196

33.7

n_cyclist_casualty_claims

1,027.0

1.8

3.2

0.0

33.0

0

253

43.5

n_pedestrian_claims

896.0

1.5

3.4

0.0

49.0

0

282

48.5

n_pedestrian_casualty_claims

715.0

1.2

2.5

0.0

37.0

0

304

52.3

population

397,237.0

683.7

432.5

0.0

4,739.0

0

3

0.5

total_roads_km

3,261.2

5.6

5.6

0.0

48.0

0

2

0.3

roads_prop_highway_arterial

67.1

0.1

0.1

0.0

0.7

2

238

41.0

no_highschool_prevalance

66.7

0.1

0.1

0.0

0.4

30

0

0.0

unemployment_rate

3,099.4

5.6

3.9

0.0

30.0

30

78

13.4

hh_avg_income

41,608,895.0

77,340.0

26,930.8

26,088.0

233,035.0

43

0

0.0

participation_rate

35,068.6

63.6

11.1

18.6

92.7

30

0

0.0

university_degree_prevalance

329.1

0.6

0.1

0.3

0.9

30

0

0.0

lone_parent_fam_prevalence

84.0

0.2

0.1

0.0

0.6

32

5

0.9

home_owner_prevalence

360.4

0.7

0.2

0.0

1.0

30

3

0.5

vandix

-130.2

-0.2

0.5

-1.5

3.0

45

0

0.0

ale_index

335.8

0.6

1.8

-2.1

7.1

28

2

0.3

canbics_index

1,170.3

2.1

1.7

0.0

8.2

28

69

11.9

Maps

claims <- ggplot() + 
  geom_sf(data = victoria_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = victoria_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = victoria_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

victoria_sp <- as(victoria_sf, "Spatial")
victoria_sp$ui <- 1:nrow(victoria_sp@data)

coords <- coordinates(victoria_sp)
victoria_nb <- poly2nb(victoria_sp, queen = TRUE)

victoria_nb
## Neighbour list object:
## Number of regions: 581 
## Number of nonzero links: 3468 
## Percentage nonzero weights: 1.02737 
## Average number of links: 5.969019
#assign nearest neighbour for no links

# victoria_nb <- assign_nearest_neighbors(victoria_nb,victoria_sp)

plot(victoria_sp, border = grey(0.5))
plot(victoria_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(victoria_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(victoria_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  victoria_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 16.303, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3913583171     -0.0017241379      0.0005813523
cyc_cc_mi <- moran.test(victoria_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  victoria_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 13.667, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3259641289     -0.0017241379      0.0005749071
pd_cc_mi <- moran.test(victoria_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  victoria_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 15.441, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3531561826     -0.0017241379      0.0005282317

BYM2 Models

nb2INLA("victoria.adj", victoria_nb)
g <- inla.read.graph(filename = "victoria.adj")


# Model Set 1: Total Casualty Claim Crashes

victoria_models_1 <- cma_models(victoria_sp@data, "Victoria", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
victoria_all <- clean_fixed_effects(victoria_models_1$fixed_effects)

map(victoria_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.2, Running = 0.866, Post = 0.148, Total = 15.3 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 2.339 0.035      2.270    2.339      2.409 2.339   0
## vandix_z    0.196 0.065      0.068    0.196      0.325 0.196   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.490 0.053      0.390    0.488      0.599 0.486
## Phi for ui       0.857 0.059      0.721    0.865      0.947 0.882
## 
## Deviance Information Criterion (DIC) ...............: 3477.53
## Deviance Information Criterion (DIC, saturated) ....: 1131.62
## Effective number of parameters .....................: 510.12
## 
## Watanabe-Akaike information criterion (WAIC) ...: 3449.37
## Effective number of parameters .................: 336.48
## 
## Marginal log-Likelihood:  -1963.43 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.2, Running = 0.894, Post = 0.162, Total = 15.2 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   2.361 0.025      2.311    2.361       2.41 2.361   0
## vandix_z      0.326 0.053      0.222    0.326       0.43 0.326   0
## ln_roads_km_c 1.237 0.073      1.095    1.237       1.38 1.237   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.642 0.050      0.547    0.640      0.744 0.639
## Phi for ui       0.981 0.019      0.930    0.987      0.999 0.996
## 
## Deviance Information Criterion (DIC) ...............: 3406.34
## Deviance Information Criterion (DIC, saturated) ....: 1060.42
## Effective number of parameters .....................: 465.82
## 
## Watanabe-Akaike information criterion (WAIC) ...: 3339.80
## Effective number of parameters .................: 285.81
## 
## Marginal log-Likelihood:  -1849.39 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.1, Running = 0.947, Post = 0.163, Total = 15.2 
## Fixed effects:
##                                mean    sd 0.025quant 0.5quant 0.975quant  mode
## (Intercept)                   2.514 0.046      2.424    2.515      2.604 2.515
## vandix_z                      0.202 0.050      0.104    0.202      0.301 0.202
## ln_roads_km_c                 1.007 0.080      0.850    1.007      1.165 1.007
## roads_prop_highway_arterial_z 0.462 0.049      0.367    0.462      0.558 0.462
## ale_index_z                   0.409 0.168      0.081    0.409      0.739 0.409
## canbics_index_z               0.340 0.156      0.032    0.340      0.646 0.340
## population_100_c              0.015 0.009     -0.002    0.015      0.031 0.015
##                               kld
## (Intercept)                     0
## vandix_z                        0
## ln_roads_km_c                   0
## roads_prop_highway_arterial_z   0
## ale_index_z                     0
## canbics_index_z                 0
## population_100_c                0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.845 0.087      0.679    0.844      1.021 0.847
## Phi for ui       0.948 0.042      0.838    0.959      0.994 0.982
## 
## Deviance Information Criterion (DIC) ...............: 3403.57
## Deviance Information Criterion (DIC, saturated) ....: 1057.64
## Effective number of parameters .....................: 453.65
## 
## Watanabe-Akaike information criterion (WAIC) ...: 3344.18
## Effective number of parameters .................: 283.52
## 
## Marginal log-Likelihood:  -1822.76 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

victoria_models_2 <- cma_models(victoria_sp@data,"Victoria","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
victoria_cyclist <- clean_fixed_effects(victoria_models_2$fixed_effects)

map(victoria_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 13.9, Running = 0.818, Post = 0.145, Total = 14.8 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) -0.141 0.070     -0.282    -0.14     -0.006 -0.14   0
## vandix_z     0.150 0.085     -0.017     0.15      0.318  0.15   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.679 0.094      0.511    0.673      0.882 0.661
## Phi for ui       0.623 0.101      0.414    0.628      0.806 0.637
## 
## Deviance Information Criterion (DIC) ...............: 1756.63
## Deviance Information Criterion (DIC, saturated) ....: 872.31
## Effective number of parameters .....................: 289.24
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1758.23
## Effective number of parameters .................: 211.24
## 
## Marginal log-Likelihood:  -721.47 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14, Running = 0.803, Post = 0.157, Total = 15 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -0.144 0.065     -0.273   -0.143     -0.018 -0.143   0
## vandix_z       0.250 0.080      0.094    0.250      0.406  0.250   0
## ln_roads_km_c  1.096 0.115      0.871    1.095      1.322  1.095   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.663 0.096      0.489    0.658      0.868 0.650
## Phi for ui       0.902 0.064      0.739    0.917      0.984 0.951
## 
## Deviance Information Criterion (DIC) ...............: 1693.39
## Deviance Information Criterion (DIC, saturated) ....: 809.07
## Effective number of parameters .....................: 240.45
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1684.94
## Effective number of parameters .................: 174.53
## 
## Marginal log-Likelihood:  -684.76 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.2, Running = 0.945, Post = 0.152, Total = 15.3 
## Fixed effects:
##                                mean    sd 0.025quant 0.5quant 0.975quant  mode
## (Intercept)                   0.039 0.080     -0.120    0.040      0.193 0.040
## vandix_z                      0.118 0.080     -0.039    0.118      0.275 0.118
## ln_roads_km_c                 0.887 0.129      0.635    0.887      1.140 0.887
## roads_prop_highway_arterial_z 0.308 0.073      0.165    0.308      0.451 0.308
## ale_index_z                   0.401 0.222     -0.033    0.401      0.838 0.401
## canbics_index_z               0.587 0.213      0.169    0.588      1.006 0.588
## population_100_c              0.031 0.014      0.004    0.031      0.058 0.031
##                               kld
## (Intercept)                     0
## vandix_z                        0
## ln_roads_km_c                   0
## roads_prop_highway_arterial_z   0
## ale_index_z                     0
## canbics_index_z                 0
## population_100_c                0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.867 0.138      0.622    0.858      1.164 0.842
## Phi for ui       0.828 0.089      0.616    0.843      0.957 0.880
## 
## Deviance Information Criterion (DIC) ...............: 1686.04
## Deviance Information Criterion (DIC, saturated) ....: 801.72
## Effective number of parameters .....................: 225.96
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1679.77
## Effective number of parameters .................: 167.01
## 
## Marginal log-Likelihood:  -686.41 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

victoria_models_3 <- cma_models(victoria_sp@data,"Victoria","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
victoria_pedestrian <- clean_fixed_effects(victoria_models_3$fixed_effects)

map(victoria_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.1, Running = 0.812, Post = 0.145, Total = 15.1 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.439 0.074     -0.589   -0.438     -0.296 -0.439   0
## vandix_z     0.350 0.084      0.184    0.349      0.516  0.349   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.633 0.099      0.458    0.626      0.847 0.615
## Phi for ui       0.791 0.096      0.569    0.805      0.937 0.837
## 
## Deviance Information Criterion (DIC) ...............: 1480.09
## Deviance Information Criterion (DIC, saturated) ....: 770.60
## Effective number of parameters .....................: 228.04
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1476.73
## Effective number of parameters .................: 167.60
## 
## Marginal log-Likelihood:  -560.34 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.2, Running = 0.812, Post = 0.15, Total = 15.2 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -0.430 0.071     -0.573   -0.429     -0.292 -0.429   0
## vandix_z       0.452 0.081      0.292    0.452      0.612  0.452   0
## ln_roads_km_c  0.854 0.118      0.625    0.854      1.086  0.854   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.619 0.091      0.458    0.613      0.816 0.602
## Phi for ui       0.950 0.047      0.822    0.964      0.996 0.988
## 
## Deviance Information Criterion (DIC) ...............: 1438.88
## Deviance Information Criterion (DIC, saturated) ....: 729.38
## Effective number of parameters .....................: 198.78
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1428.78
## Effective number of parameters .................: 144.38
## 
## Marginal log-Likelihood:  -540.86 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.2, Running = 1.56, Post = 1.14, Total = 16.9 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -0.313 0.090     -0.493   -0.312     -0.141
## vandix_z                       0.317 0.082      0.156    0.317      0.478
## ln_roads_km_c                  0.578 0.135      0.313    0.578      0.845
## roads_prop_highway_arterial_z  0.384 0.079      0.228    0.383      0.539
## ale_index_z                    0.540 0.237      0.080    0.539      1.008
## canbics_index_z                0.326 0.235     -0.136    0.327      0.786
## population_100_c               0.037 0.014      0.010    0.037      0.064
##                                 mode kld
## (Intercept)                   -0.312   0
## vandix_z                       0.317   0
## ln_roads_km_c                  0.578   0
## roads_prop_highway_arterial_z  0.383   0
## ale_index_z                    0.539   0
## canbics_index_z                0.327   0
## population_100_c               0.037   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.758 0.122      0.543    0.750      1.023 0.734
## Phi for ui       0.924 0.066      0.745    0.944      0.993 0.980
## 
## Deviance Information Criterion (DIC) ...............: 1426.26
## Deviance Information Criterion (DIC, saturated) ....: 716.76
## Effective number of parameters .....................: 185.00
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1417.96
## Effective number of parameters .................: 136.61
## 
## Marginal log-Likelihood:  -542.77 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
victoria_results <- bind_rows(victoria_all %>% mutate(Region = "Victoria",Outcome = "All Injury Claims"),
                         victoria_cyclist %>% mutate(Region = "Victoria",Outcome = "Cyclist Injury Claims"),
                         victoria_pedestrian %>% mutate(Region = "Victoria",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

Central Island - Powell River

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Central Island-Powell River

n_claims

611

62,992.0

103.1

164.7

0.0

1,934.0

0

17

2.8

n_casualty_claims

10,762.0

17.6

28.3

0.0

203.0

0

47

7.7

n_cyclist_claims

475.0

0.8

1.5

0.0

14.0

0

365

59.7

n_cyclist_casualty_claims

271.0

0.4

0.9

0.0

10.0

0

435

71.2

n_pedestrian_claims

673.0

1.1

2.4

0.0

23.0

0

361

59.1

n_pedestrian_casualty_claims

547.0

0.9

2.0

0.0

18.0

0

387

63.3

population

357,193.0

584.6

307.8

0.0

2,466.0

0

14

2.3

total_roads_km

5,652.1

9.3

9.9

0.0

86.5

0

10

1.6

roads_prop_highway_arterial

78.2

0.1

0.2

0.0

0.8

10

262

42.9

no_highschool_prevalance

101.2

0.2

0.1

0.0

0.6

38

1

0.2

unemployment_rate

4,784.2

8.3

5.2

0.0

40.0

38

43

7.0

hh_avg_income

36,682,051.0

65,738.4

17,197.1

29,925.0

211,531.0

53

0

0.0

participation_rate

32,393.0

56.5

11.0

12.0

86.4

38

0

0.0

university_degree_prevalance

300.2

0.5

0.1

0.2

0.8

38

0

0.0

lone_parent_fam_prevalence

95.7

0.2

0.1

0.0

0.6

38

3

0.5

home_owner_prevalence

426.2

0.7

0.2

0.0

1.0

38

1

0.2

vandix

111.8

0.2

0.7

-1.1

4.0

53

0

0.0

ale_index

-590.4

-1.1

0.7

-2.1

1.5

57

0

0.0

canbics_index

326.8

0.6

0.9

0.0

4.2

55

241

39.4

Maps

claims <- ggplot() + 
  geom_sf(data = nanaimo_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = nanaimo_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = nanaimo_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

nanaimo_sp <- as(nanaimo_sf, "Spatial")
nanaimo_sp$ui <- 1:nrow(nanaimo_sp@data)

coords <- coordinates(nanaimo_sp)
nanaimo_nb <- poly2nb(nanaimo_sp, queen = TRUE)
## Warning in poly2nb(nanaimo_sp, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in poly2nb(nanaimo_sp, queen = TRUE): neighbour object has 6 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
nanaimo_nb
## Neighbour list object:
## Number of regions: 611 
## Number of nonzero links: 3330 
## Percentage nonzero weights: 0.8919938 
## Average number of links: 5.450082 
## 1 region with no links:
## 169
## 6 disjoint connected subgraphs
#assign nearest neighbour for no links

nanaimo_nb <- assign_nearest_neighbors(nanaimo_nb,nanaimo_sp)

plot(nanaimo_sp, border = grey(0.5))
plot(nanaimo_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(nanaimo_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(nanaimo_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  nanaimo_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 12.06, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3077074730     -0.0016393443      0.0006579567
cyc_cc_mi <- moran.test(nanaimo_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  nanaimo_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 6.5732, p-value = 2.463e-11
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1652016950     -0.0016393443      0.0006442558
pd_cc_mi <- moran.test(nanaimo_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  nanaimo_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 8.4624, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2136910003     -0.0016393443      0.0006474725

BYM2 Models

nb2INLA("nanaimo.adj", nanaimo_nb)
g <- inla.read.graph(filename = "nanaimo.adj")


# Model Set 1: Total Casualty Claim Crashes

nanaimo_models_1 <- cma_models(nanaimo_sp@data, "Central Island-Powell River", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
nanaimo_all <- clean_fixed_effects(nanaimo_models_1$fixed_effects)

map(nanaimo_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 15.4, Running = 0.911, Post = 0.19, Total = 16.5 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 1.926 0.050      1.826    1.926      2.024 1.926   0
## vandix_z    0.295 0.061      0.175    0.295      0.416 0.295   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.410 0.058      0.309    0.406      0.537 0.397
## Phi for ui       0.571 0.092      0.383    0.574      0.741 0.583
## 
## Deviance Information Criterion (DIC) ...............: 3536.92
## Deviance Information Criterion (DIC, saturated) ....: 1244.46
## Effective number of parameters .....................: 555.03
## 
## Watanabe-Akaike information criterion (WAIC) ...: 3574.86
## Effective number of parameters .................: 405.59
## 
## Marginal log-Likelihood:  -2193.01 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 15.9, Running = 1.18, Post = 0.168, Total = 17.2 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   1.357 0.060      1.238    1.357      1.475 1.357   0
## vandix_z      0.435 0.053      0.330    0.435      0.539 0.435   0
## ln_roads_km_c 1.024 0.069      0.888    1.023      1.160 1.023   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.442 0.054      0.345    0.438      0.558 0.432
## Phi for ui       0.736 0.057      0.613    0.740      0.836 0.748
## 
## Deviance Information Criterion (DIC) ...............: 3467.76
## Deviance Information Criterion (DIC, saturated) ....: 1175.29
## Effective number of parameters .....................: 523.49
## 
## Watanabe-Akaike information criterion (WAIC) ...: 3456.68
## Effective number of parameters .................: 356.15
## 
## Marginal log-Likelihood:  -2098.24 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 15.6, Running = 1.08, Post = 0.193, Total = 16.9 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                    1.961 0.117      1.730    1.962      2.190
## vandix_z                       0.317 0.046      0.227    0.317      0.408
## ln_roads_km_c                  0.650 0.073      0.506    0.650      0.794
## roads_prop_highway_arterial_z  0.668 0.047      0.575    0.668      0.761
## ale_index_z                    1.419 0.246      0.937    1.419      1.903
## canbics_index_z               -0.780 0.231     -1.233   -0.779     -0.328
## population_100_c               0.063 0.013      0.037    0.063      0.090
##                                 mode kld
## (Intercept)                    1.962   0
## vandix_z                       0.317   0
## ln_roads_km_c                  0.650   0
## roads_prop_highway_arterial_z  0.668   0
## ale_index_z                    1.419   0
## canbics_index_z               -0.779   0
## population_100_c               0.063   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.670 0.085      0.517    0.665      0.852 0.656
## Phi for ui       0.726 0.060      0.598    0.730      0.832 0.738
## 
## Deviance Information Criterion (DIC) ...............: 3436.59
## Deviance Information Criterion (DIC, saturated) ....: 1144.12
## Effective number of parameters .....................: 492.67
## 
## Watanabe-Akaike information criterion (WAIC) ...: 3402.21
## Effective number of parameters .................: 324.32
## 
## Marginal log-Likelihood:  -2024.42 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

nanaimo_models_2 <- cma_models(nanaimo_sp@data,"Central Island-Powell River","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
nanaimo_cyclist <- clean_fixed_effects(nanaimo_models_2$fixed_effects)

map(nanaimo_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 15.6, Running = 0.969, Post = 0.154, Total = 16.8 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.429 0.110     -1.655   -1.427     -1.220 -1.427   0
## vandix_z     0.335 0.075      0.188    0.335      0.483  0.335   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                  mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.06 0.236      0.680    1.031      1.605 0.972
## Phi for ui       0.32 0.142      0.089    0.305      0.626 0.266
## 
## Deviance Information Criterion (DIC) ...............: 1012.10
## Deviance Information Criterion (DIC, saturated) ....: 614.45
## Effective number of parameters .....................: 137.56
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1008.51
## Effective number of parameters .................: 106.42
## 
## Marginal log-Likelihood:  -435.69 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 15.3, Running = 1.01, Post = 0.158, Total = 16.5 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -1.731 0.139     -2.012   -1.729     -1.465 -1.729   0
## vandix_z       0.399 0.076      0.250    0.399      0.548  0.399   0
## ln_roads_km_c  0.506 0.117      0.280    0.505      0.737  0.505   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.927 0.231      0.561    0.898      1.463 0.841
## Phi for ui       0.551 0.151      0.249    0.557      0.823 0.577
## 
## Deviance Information Criterion (DIC) ...............: 993.69
## Deviance Information Criterion (DIC, saturated) ....: 596.04
## Effective number of parameters .....................: 126.27
## 
## Watanabe-Akaike information criterion (WAIC) ...: 993.42
## Effective number of parameters .................: 100.60
## 
## Marginal log-Likelihood:  -431.79 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 15.2, Running = 1.09, Post = 0.173, Total = 16.5 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -1.152 0.208     -1.568   -1.149     -0.753
## vandix_z                       0.354 0.077      0.203    0.354      0.504
## ln_roads_km_c                  0.292 0.140      0.018    0.291      0.568
## roads_prop_highway_arterial_z  0.431 0.082      0.272    0.431      0.592
## ale_index_z                    1.156 0.413      0.350    1.154      1.969
## canbics_index_z               -0.321 0.366     -1.042   -0.319      0.394
## population_100_c               0.060 0.024      0.014    0.060      0.107
##                                 mode kld
## (Intercept)                   -1.149   0
## vandix_z                       0.354   0
## ln_roads_km_c                  0.291   0
## roads_prop_highway_arterial_z  0.431   0
## ale_index_z                    1.154   0
## canbics_index_z               -0.319   0
## population_100_c               0.060   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision for ui 1.177 0.313      0.691    1.133      1.915 1.05
## Phi for ui       0.496 0.156      0.197    0.497      0.788 0.51
## 
## Deviance Information Criterion (DIC) ...............: 971.41
## Deviance Information Criterion (DIC, saturated) ....: 573.76
## Effective number of parameters .....................: 113.46
## 
## Watanabe-Akaike information criterion (WAIC) ...: 972.04
## Effective number of parameters .................: 92.21
## 
## Marginal log-Likelihood:  -433.56 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

nanaimo_models_3 <- cma_models(nanaimo_sp@data,"Central Island-Powell River","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
nanaimo_pedestrian <- clean_fixed_effects(nanaimo_models_3$fixed_effects)

map(nanaimo_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 15.2, Running = 1.02, Post = 0.159, Total = 16.3 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.252 0.116     -1.494   -1.248     -1.036 -1.248   0
## vandix_z     0.556 0.078      0.405    0.556      0.710  0.556   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.605 0.100      0.439    0.595      0.830 0.572
## Phi for ui       0.206 0.128      0.030    0.180      0.511 0.102
## 
## Deviance Information Criterion (DIC) ...............: 1356.86
## Deviance Information Criterion (DIC, saturated) ....: 791.04
## Effective number of parameters .....................: 274.96
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1378.25
## Effective number of parameters .................: 208.18
## 
## Marginal log-Likelihood:  -648.63 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 15.5, Running = 1.05, Post = 0.162, Total = 16.8 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -1.598 0.140     -1.882   -1.595     -1.330 -1.595   0
## vandix_z       0.611 0.080      0.456    0.611      0.769  0.611   0
## ln_roads_km_c  0.653 0.129      0.401    0.652      0.908  0.652   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.468 0.102      0.302    0.456      0.701 0.433
## Phi for ui       0.587 0.140      0.295    0.596      0.830 0.627
## 
## Deviance Information Criterion (DIC) ...............: 1314.98
## Deviance Information Criterion (DIC, saturated) ....: 749.16
## Effective number of parameters .....................: 241.00
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1327.06
## Effective number of parameters .................: 182.32
## 
## Marginal log-Likelihood:  -641.60 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 15.4, Running = 1.04, Post = 0.155, Total = 16.5 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -0.755 0.205     -1.167   -0.752     -0.361
## vandix_z                       0.518 0.075      0.370    0.518      0.666
## ln_roads_km_c                  0.328 0.141      0.052    0.327      0.605
## roads_prop_highway_arterial_z  0.552 0.081      0.394    0.551      0.713
## ale_index_z                    2.371 0.393      1.605    2.369      3.147
## canbics_index_z               -1.147 0.349     -1.838   -1.146     -0.468
## population_100_c               0.115 0.023      0.071    0.115      0.160
##                                 mode kld
## (Intercept)                   -0.752   0
## vandix_z                       0.518   0
## ln_roads_km_c                  0.327   0
## roads_prop_highway_arterial_z  0.551   0
## ale_index_z                    2.369   0
## canbics_index_z               -1.145   0
## population_100_c               0.115   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.749 0.161      0.489    0.731      1.119 0.692
## Phi for ui       0.446 0.150      0.166    0.444      0.736 0.451
## 
## Deviance Information Criterion (DIC) ...............: 1265.20
## Deviance Information Criterion (DIC, saturated) ....: 699.38
## Effective number of parameters .....................: 201.41
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1265.47
## Effective number of parameters .................: 150.79
## 
## Marginal log-Likelihood:  -617.89 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
nanaimo_results <- bind_rows(nanaimo_all %>% mutate(Region = "Central Island-Powell River",Outcome = "All Injury Claims"),
                         nanaimo_cyclist %>% mutate(Region = "Central Island-Powell River",Outcome = "Cyclist Injury Claims"),
                         nanaimo_pedestrian %>% mutate(Region = "Central Island-Powell River",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

Okanagan

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Okanagan

n_claims

466

70,024.0

150.3

325.5

0.0

4,324.0

0

10

2.1

n_casualty_claims

13,236.0

28.4

66.6

0.0

938.0

0

38

8.2

n_cyclist_claims

761.0

1.6

3.8

0.0

38.0

0

252

54.1

n_cyclist_casualty_claims

485.0

1.0

2.3

0.0

23.0

0

286

61.4

n_pedestrian_claims

696.0

1.5

3.9

0.0

41.0

0

259

55.6

n_pedestrian_casualty_claims

564.0

1.2

3.2

0.0

35.0

0

274

58.8

population

336,628.0

722.4

489.5

0.0

4,622.0

0

10

2.1

total_roads_km

4,759.1

10.2

13.6

0.0

139.3

0

2

0.4

roads_prop_highway_arterial

60.1

0.1

0.2

0.0

1.0

2

199

42.7

no_highschool_prevalance

62.0

0.2

0.1

0.0

0.4

88

4

0.9

unemployment_rate

2,933.5

7.8

4.7

0.0

33.3

88

28

6.0

hh_avg_income

26,784,003.0

72,980.9

26,808.0

26,241.0

247,493.0

99

0

0.0

participation_rate

23,166.9

61.3

12.9

9.2

91.5

88

0

0.0

university_degree_prevalance

199.0

0.5

0.1

0.3

0.8

88

0

0.0

lone_parent_fam_prevalence

59.9

0.2

0.1

0.0

0.5

88

4

0.9

home_owner_prevalence

275.1

0.7

0.2

0.0

1.3

88

0

0.0

vandix

42.6

0.1

0.6

-1.2

2.2

99

0

0.0

ale_index

-314.5

-0.8

1.0

-2.1

2.3

93

1

0.2

canbics_index

709.8

1.9

1.8

0.0

7.2

93

82

17.6

Maps

claims <- ggplot() + 
  geom_sf(data = okanagan_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = okanagan_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = okanagan_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

okanagan_sp <- as(okanagan_sf, "Spatial")
okanagan_sp$ui <- 1:nrow(okanagan_sp@data)

coords <- coordinates(okanagan_sp)
okanagan_nb <- poly2nb(okanagan_sp, queen = TRUE,snap = 1)

okanagan_nb
## Neighbour list object:
## Number of regions: 466 
## Number of nonzero links: 2728 
## Percentage nonzero weights: 1.25624 
## Average number of links: 5.854077
#assign nearest neighbour for no links

# okanagan_nb <- assign_nearest_neighbors(okanagan_nb,okanagan_sp)

plot(okanagan_sp, border = grey(0.5))
plot(okanagan_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(okanagan_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(okanagan_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  okanagan_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 14.429, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3626842832     -0.0021505376      0.0006393647
cyc_cc_mi <- moran.test(okanagan_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  okanagan_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 14.984, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.4004592704     -0.0021505376      0.0007220019
pd_cc_mi <- moran.test(okanagan_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  okanagan_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 8.8431, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2316536313     -0.0021505376      0.0006990321

BYM2 Models

nb2INLA("okanagan.adj", okanagan_nb)
g <- inla.read.graph(filename = "okanagan.adj")


# Model Set 1: Total Casualty Claim Crashes

okanagan_models_1 <- cma_models(okanagan_sp@data, "Okanagan", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
okanagan_all <- clean_fixed_effects(okanagan_models_1$fixed_effects)

map(okanagan_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.1, Running = 0.829, Post = 0.166, Total = 15.1 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 2.164 0.047      2.072    2.164      2.255 2.164   0
## vandix_z    0.344 0.092      0.163    0.344      0.525 0.344   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.237 0.037      0.172    0.234      0.317 0.230
## Phi for ui       0.843 0.059      0.704    0.851      0.934 0.868
## 
## Deviance Information Criterion (DIC) ...............: 2782.89
## Deviance Information Criterion (DIC, saturated) ....: 931.94
## Effective number of parameters .....................: 422.95
## 
## Watanabe-Akaike information criterion (WAIC) ...: 2786.87
## Effective number of parameters .................: 293.18
## 
## Marginal log-Likelihood:  -1785.42 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.3, Running = 0.919, Post = 0.134, Total = 15.4 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   1.579 0.053      1.476    1.579      1.682 1.579   0
## vandix_z      0.464 0.074      0.318    0.463      0.610 0.463   0
## ln_roads_km_c 1.145 0.073      1.002    1.145      1.289 1.145   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.241 0.022       0.20    0.241      0.286 0.241
## Phi for ui       0.985 0.013       0.95    0.989      0.998 0.996
## 
## Deviance Information Criterion (DIC) ...............: 2705.48
## Deviance Information Criterion (DIC, saturated) ....: 854.52
## Effective number of parameters .....................: 386.80
## 
## Watanabe-Akaike information criterion (WAIC) ...: 2655.01
## Effective number of parameters .................: 237.31
## 
## Marginal log-Likelihood:  -1689.24 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.2, Running = 0.983, Post = 0.153, Total = 15.3 
## Fixed effects:
##                                mean    sd 0.025quant 0.5quant 0.975quant  mode
## (Intercept)                   1.962 0.084      1.796    1.962      2.126 1.962
## vandix_z                      0.355 0.065      0.228    0.355      0.484 0.355
## ln_roads_km_c                 0.949 0.076      0.801    0.949      1.098 0.949
## roads_prop_highway_arterial_z 0.605 0.054      0.499    0.605      0.711 0.605
## ale_index_z                   0.510 0.263     -0.006    0.510      1.026 0.510
## canbics_index_z               0.016 0.215     -0.406    0.016      0.438 0.016
## population_100_c              0.031 0.009      0.014    0.031      0.048 0.031
##                               kld
## (Intercept)                     0
## vandix_z                        0
## ln_roads_km_c                   0
## roads_prop_highway_arterial_z   0
## ale_index_z                     0
## canbics_index_z                 0
## population_100_c                0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.339 0.036      0.271    0.339      0.412 0.340
## Phi for ui       0.979 0.019      0.929    0.985      0.998 0.995
## 
## Deviance Information Criterion (DIC) ...............: 2695.73
## Deviance Information Criterion (DIC, saturated) ....: 844.77
## Effective number of parameters .....................: 371.48
## 
## Watanabe-Akaike information criterion (WAIC) ...: 2644.06
## Effective number of parameters .................: 228.01
## 
## Marginal log-Likelihood:  -1650.75 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

okanagan_models_2 <- cma_models(okanagan_sp@data,"Okanagan","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
okanagan_cyclist <- clean_fixed_effects(okanagan_models_2$fixed_effects)

map(okanagan_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.1, Running = 0.748, Post = 0.142, Total = 15 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.113 0.110     -1.337   -1.111     -0.903 -1.112   0
## vandix_z     0.642 0.109      0.429    0.642      0.858  0.642   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.426 0.093      0.269    0.417      0.634 0.402
## Phi for ui       0.788 0.095      0.568    0.802      0.933 0.834
## 
## Deviance Information Criterion (DIC) ...............: 995.71
## Deviance Information Criterion (DIC, saturated) ....: 526.98
## Effective number of parameters .....................: 160.92
## 
## Watanabe-Akaike information criterion (WAIC) ...: 993.54
## Effective number of parameters .................: 117.40
## 
## Marginal log-Likelihood:  -488.26 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.1, Running = 0.832, Post = 0.155, Total = 15.1 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -1.507 0.140     -1.787   -1.505     -1.239 -1.505   0
## vandix_z       0.705 0.105      0.500    0.704      0.913  0.704   0
## ln_roads_km_c  0.721 0.121      0.484    0.721      0.961  0.721   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.339 0.062       0.23    0.334      0.474 0.327
## Phi for ui       0.951 0.041       0.84    0.963      0.995 0.985
## 
## Deviance Information Criterion (DIC) ...............: 964.03
## Deviance Information Criterion (DIC, saturated) ....: 495.30
## Effective number of parameters .....................: 140.52
## 
## Watanabe-Akaike information criterion (WAIC) ...: 957.29
## Effective number of parameters .................: 101.14
## 
## Marginal log-Likelihood:  -477.10 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.1, Running = 0.86, Post = 0.162, Total = 15.1 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -1.014 0.165     -1.346   -1.012     -0.697
## vandix_z                       0.566 0.101      0.369    0.566      0.765
## ln_roads_km_c                  0.616 0.137      0.348    0.616      0.887
## roads_prop_highway_arterial_z  0.395 0.089      0.221    0.394      0.569
## ale_index_z                    0.703 0.379     -0.040    0.703      1.448
## canbics_index_z                0.354 0.292     -0.220    0.355      0.925
## population_100_c               0.036 0.015      0.008    0.036      0.065
##                                 mode kld
## (Intercept)                   -1.012   0
## vandix_z                       0.566   0
## ln_roads_km_c                  0.616   0
## roads_prop_highway_arterial_z  0.394   0
## ale_index_z                    0.703   0
## canbics_index_z                0.355   0
## population_100_c               0.036   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.428 0.082      0.288    0.421      0.608 0.409
## Phi for ui       0.961 0.038      0.858    0.973      0.997 0.991
## 
## Deviance Information Criterion (DIC) ...............: 951.76
## Deviance Information Criterion (DIC, saturated) ....: 483.03
## Effective number of parameters .....................: 124.22
## 
## Watanabe-Akaike information criterion (WAIC) ...: 945.65
## Effective number of parameters .................: 91.39
## 
## Marginal log-Likelihood:  -480.27 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

okanagan_models_3 <- cma_models(okanagan_sp@data,"Okanagan","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
okanagan_pedestrian <- clean_fixed_effects(okanagan_models_3$fixed_effects)

map(okanagan_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.2, Running = 0.777, Post = 0.143, Total = 15.1 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.042 0.117     -1.284   -1.038     -0.822 -1.039   0
## vandix_z     0.668 0.110      0.452    0.668      0.884  0.668   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.525 0.088      0.374    0.517      0.720 0.501
## Phi for ui       0.328 0.123      0.121    0.317      0.593 0.291
## 
## Deviance Information Criterion (DIC) ...............: 1114.17
## Deviance Information Criterion (DIC, saturated) ....: 620.03
## Effective number of parameters .....................: 220.56
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1132.39
## Effective number of parameters .................: 167.12
## 
## Marginal log-Likelihood:  -554.86 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.2, Running = 0.809, Post = 0.157, Total = 15.2 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -1.439 0.138     -1.716   -1.437     -1.173 -1.437   0
## vandix_z       0.731 0.119      0.497    0.730      0.965  0.730   0
## ln_roads_km_c  0.805 0.135      0.541    0.805      1.073  0.805   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.372 0.084      0.232    0.364      0.561 0.349
## Phi for ui       0.751 0.114      0.489    0.768      0.925 0.808
## 
## Deviance Information Criterion (DIC) ...............: 1070.51
## Deviance Information Criterion (DIC, saturated) ....: 576.37
## Effective number of parameters .....................: 185.13
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1071.08
## Effective number of parameters .................: 135.15
## 
## Marginal log-Likelihood:  -542.90 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.4, Running = 0.914, Post = 0.149, Total = 15.4 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -0.861 0.160     -1.182   -0.859     -0.554
## vandix_z                       0.558 0.113      0.338    0.558      0.780
## ln_roads_km_c                  0.632 0.150      0.337    0.632      0.927
## roads_prop_highway_arterial_z  0.575 0.088      0.404    0.575      0.748
## ale_index_z                    1.294 0.391      0.531    1.293      2.065
## canbics_index_z               -0.245 0.321     -0.877   -0.244      0.384
## population_100_c               0.049 0.016      0.019    0.049      0.081
##                                 mode kld
## (Intercept)                   -0.859   0
## vandix_z                       0.558   0
## ln_roads_km_c                  0.632   0
## roads_prop_highway_arterial_z  0.575   0
## ale_index_z                    1.293   0
## canbics_index_z               -0.244   0
## population_100_c               0.049   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.508 0.122      0.310    0.494      0.786 0.468
## Phi for ui       0.707 0.132      0.408    0.725      0.911 0.772
## 
## Deviance Information Criterion (DIC) ...............: 1040.50
## Deviance Information Criterion (DIC, saturated) ....: 546.36
## Effective number of parameters .....................: 164.38
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1035.97
## Effective number of parameters .................: 118.86
## 
## Marginal log-Likelihood:  -534.59 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
okanagan_results <- bind_rows(okanagan_all %>% mutate(Region = "Okanagan",Outcome = "All Injury Claims"),
                         okanagan_cyclist %>% mutate(Region = "Okanagan",Outcome = "Cyclist Injury Claims"),
                         okanagan_pedestrian %>% mutate(Region = "Okanagan",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

Fraser Valley

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Fraser Valley

n_claims

472

68,269.0

144.6

279.9

0.0

2,538.0

0

4

0.8

n_casualty_claims

15,730.0

33.3

65.2

0.0

556.0

0

17

3.6

n_cyclist_claims

581.0

1.2

2.4

0.0

21.0

0

248

52.5

n_cyclist_casualty_claims

360.0

0.8

1.6

0.0

12.0

0

302

64.0

n_pedestrian_claims

894.0

1.9

3.7

0.0

32.0

0

205

43.4

n_pedestrian_casualty_claims

725.0

1.5

3.0

0.0

25.0

0

225

47.7

population

309,493.0

655.7

341.8

0.0

2,550.0

0

4

0.8

total_roads_km

3,678.9

7.8

9.6

0.1

81.8

0

0

0.0

roads_prop_highway_arterial

40.5

0.1

0.1

0.0

0.7

0

281

59.5

no_highschool_prevalance

82.9

0.2

0.1

0.1

0.5

88

0

0.0

unemployment_rate

2,706.8

7.0

5.1

0.0

30.0

88

45

9.5

hh_avg_income

27,748,698.0

75,404.1

21,091.5

24,502.0

161,236.0

104

0

0.0

participation_rate

24,427.3

63.6

11.6

12.7

82.2

88

0

0.0

university_degree_prevalance

170.5

0.4

0.1

0.1

0.7

88

0

0.0

lone_parent_fam_prevalence

63.8

0.2

0.1

0.0

1.0

87

3

0.6

home_owner_prevalence

284.5

0.7

0.2

0.0

1.0

88

1

0.2

vandix

112.4

0.3

0.6

-0.9

2.0

104

0

0.0

ale_index

-89.0

-0.2

1.7

-2.1

6.7

101

1

0.2

canbics_index

514.5

1.4

1.3

0.0

4.9

101

81

17.2

Maps

claims <- ggplot() + 
  geom_sf(data = fraser_valley_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = fraser_valley_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = fraser_valley_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

fraser_valley_sp <- as(fraser_valley_sf, "Spatial")
fraser_valley_sp$ui <- 1:nrow(fraser_valley_sp@data)

coords <- coordinates(fraser_valley_sp)
fraser_valley_nb <- poly2nb(fraser_valley_sp, queen = TRUE)
## Warning in poly2nb(fraser_valley_sp, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
fraser_valley_nb
## Neighbour list object:
## Number of regions: 472 
## Number of nonzero links: 2822 
## Percentage nonzero weights: 1.266698 
## Average number of links: 5.978814 
## 2 disjoint connected subgraphs
#assign nearest neighbour for no links

fraser_valley_nb <- assign_nearest_neighbors(fraser_valley_nb,fraser_valley_sp)

plot(fraser_valley_sp, border = grey(0.5))
plot(fraser_valley_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(fraser_valley_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(fraser_valley_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  fraser_valley_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 8.2045, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2194714130     -0.0021231423      0.0007294782
cyc_cc_mi <- moran.test(fraser_valley_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  fraser_valley_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 8.2775, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2236414065     -0.0021231423      0.0007438924
pd_cc_mi <- moran.test(fraser_valley_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  fraser_valley_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 9.3134, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2508545992     -0.0021231423      0.0007378148

BYM2 Models

nb2INLA("fraser_valley.adj", fraser_valley_nb)
g <- inla.read.graph(filename = "fraser_valley.adj")


# Model Set 1: Total Casualty Claim Crashes

fraser_valley_models_1 <- cma_models(fraser_valley_sp@data, "Fraser Valley", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
fraser_valley_all <- clean_fixed_effects(fraser_valley_models_1$fixed_effects)

map(fraser_valley_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.4, Running = 0.907, Post = 0.133, Total = 15.4 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 2.479 0.040      2.399    2.479      2.557 2.479   0
## vandix_z    0.232 0.085      0.065    0.232      0.399 0.232   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.299 0.022      0.258    0.298      0.344 0.297
## Phi for ui       0.767 0.354      0.001    0.985      1.000 1.000
## 
## Deviance Information Criterion (DIC) ...............: 2987.34
## Deviance Information Criterion (DIC, saturated) ....: 924.07
## Effective number of parameters .....................: 430.45
## 
## Watanabe-Akaike information criterion (WAIC) ...: 2950.25
## Effective number of parameters .................: 274.19
## 
## Marginal log-Likelihood:  -1818.76 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.4, Running = 0.84, Post = 0.14, Total = 15.4 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   2.265 0.043      2.180    2.265      2.350 2.265   0
## vandix_z      0.285 0.066      0.156    0.285      0.415 0.285   0
## ln_roads_km_c 1.104 0.067      0.972    1.104      1.236 1.104   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.712 0.069      0.585    0.709      0.856 0.703
## Phi for ui       0.802 0.055      0.680    0.807      0.893 0.818
## 
## Deviance Information Criterion (DIC) ...............: 2985.16
## Deviance Information Criterion (DIC, saturated) ....: 921.90
## Effective number of parameters .....................: 417.80
## 
## Watanabe-Akaike information criterion (WAIC) ...: 2950.49
## Effective number of parameters .................: 269.04
## 
## Marginal log-Likelihood:  -1674.82 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.5, Running = 0.949, Post = 0.184, Total = 15.6 
## Fixed effects:
##                                mean    sd 0.025quant 0.5quant 0.975quant  mode
## (Intercept)                   2.719 0.060      2.601    2.719      2.837 2.719
## vandix_z                      0.140 0.060      0.022    0.139      0.257 0.139
## ln_roads_km_c                 0.862 0.072      0.721    0.862      1.004 0.862
## roads_prop_highway_arterial_z 0.580 0.059      0.464    0.580      0.695 0.580
## ale_index_z                   0.245 0.127     -0.004    0.245      0.494 0.245
## canbics_index_z               0.183 0.142     -0.097    0.183      0.458 0.183
## population_100_c              0.027 0.012      0.003    0.027      0.051 0.027
##                               kld
## (Intercept)                     0
## vandix_z                        0
## ln_roads_km_c                   0
## roads_prop_highway_arterial_z   0
## ale_index_z                     0
## canbics_index_z                 0
## population_100_c                0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.953 0.104      0.760    0.949      1.170 0.944
## Phi for ui       0.805 0.066      0.658    0.812      0.913 0.825
## 
## Deviance Information Criterion (DIC) ...............: 2973.51
## Deviance Information Criterion (DIC, saturated) ....: 910.24
## Effective number of parameters .....................: 403.86
## 
## Watanabe-Akaike information criterion (WAIC) ...: 2931.21
## Effective number of parameters .................: 256.46
## 
## Marginal log-Likelihood:  -1640.90 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

fraser_valley_models_2 <- cma_models(fraser_valley_sp@data,"Fraser Valley","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
fraser_valley_cyclist <- clean_fixed_effects(fraser_valley_models_2$fixed_effects)

map(fraser_valley_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.4, Running = 0.852, Post = 0.136, Total = 15.4 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.220 0.125     -1.478   -1.216     -0.984 -1.216   0
## vandix_z     0.338 0.113      0.113    0.338      0.557  0.338   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.679 0.113      0.485    0.669      0.930 0.648
## Phi for ui       0.315 0.130      0.098    0.303      0.593 0.275
## 
## Deviance Information Criterion (DIC) ...............: 997.87
## Deviance Information Criterion (DIC, saturated) ....: 582.62
## Effective number of parameters .....................: 185.63
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1003.75
## Effective number of parameters .................: 138.73
## 
## Marginal log-Likelihood:  -315.14 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.4, Running = 1.03, Post = 0.189, Total = 15.6 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -1.360 0.124     -1.610   -1.358     -1.123 -1.358   0
## vandix_z       0.286 0.114      0.063    0.286      0.509  0.286   0
## ln_roads_km_c  0.753 0.127      0.505    0.752      1.004  0.752   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.623 0.114      0.427    0.614      0.874 0.596
## Phi for ui       0.752 0.120      0.479    0.769      0.934 0.813
## 
## Deviance Information Criterion (DIC) ...............: 961.35
## Deviance Information Criterion (DIC, saturated) ....: 546.10
## Effective number of parameters .....................: 153.84
## 
## Watanabe-Akaike information criterion (WAIC) ...: 964.81
## Effective number of parameters .................: 117.69
## 
## Marginal log-Likelihood:  -303.34 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.9, Running = 1.03, Post = 0.234, Total = 16.1 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -0.791 0.145     -1.082   -0.789     -0.513
## vandix_z                       0.081 0.112     -0.139    0.081      0.299
## ln_roads_km_c                  0.467 0.133      0.207    0.466      0.731
## roads_prop_highway_arterial_z  0.601 0.101      0.404    0.600      0.800
## ale_index_z                    0.695 0.214      0.276    0.695      1.115
## canbics_index_z                0.091 0.299     -0.494    0.091      0.678
## population_100_c               0.046 0.023      0.001    0.046      0.091
##                                 mode kld
## (Intercept)                   -0.789   0
## vandix_z                       0.081   0
## ln_roads_km_c                  0.466   0
## roads_prop_highway_arterial_z  0.600   0
## ale_index_z                    0.695   0
## canbics_index_z                0.091   0
## population_100_c               0.046   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.897 0.186      0.588    0.879      1.316 0.842
## Phi for ui       0.611 0.160      0.278    0.623      0.882 0.659
## 
## Deviance Information Criterion (DIC) ...............: 941.71
## Deviance Information Criterion (DIC, saturated) ....: 526.46
## Effective number of parameters .....................: 136.45
## 
## Watanabe-Akaike information criterion (WAIC) ...: 942.40
## Effective number of parameters .................: 105.00
## 
## Marginal log-Likelihood:  -300.28 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

fraser_valley_models_3 <- cma_models(fraser_valley_sp@data,"Fraser Valley","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
fraser_valley_pedestrian <- clean_fixed_effects(fraser_valley_models_3$fixed_effects)

map(fraser_valley_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.6, Running = 0.928, Post = 0.165, Total = 15.7 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.585 0.096     -0.779   -0.583     -0.402 -0.583   0
## vandix_z     0.503 0.095      0.318    0.503      0.689  0.503   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.740 0.104      0.558    0.732      0.965 0.716
## Phi for ui       0.372 0.115      0.166    0.366      0.610 0.356
## 
## Deviance Information Criterion (DIC) ...............: 1347.02
## Deviance Information Criterion (DIC, saturated) ....: 695.41
## Effective number of parameters .....................: 236.66
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1353.70
## Effective number of parameters .................: 174.42
## 
## Marginal log-Likelihood:  -517.06 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.6, Running = 0.821, Post = 0.146, Total = 15.5 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -0.730 0.092     -0.915   -0.729     -0.552 -0.729   0
## vandix_z       0.447 0.091      0.269    0.447      0.627  0.447   0
## ln_roads_km_c  0.813 0.104      0.608    0.813      1.017  0.813   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.660 0.103      0.476    0.653      0.882 0.642
## Phi for ui       0.852 0.087      0.638    0.869      0.970 0.912
## 
## Deviance Information Criterion (DIC) ...............: 1295.95
## Deviance Information Criterion (DIC, saturated) ....: 644.34
## Effective number of parameters .....................: 194.13
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1293.89
## Effective number of parameters .................: 143.01
## 
## Marginal log-Likelihood:  -495.02 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.5, Running = 0.87, Post = 0.169, Total = 15.5 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -0.210 0.114     -0.437   -0.208      0.010
## vandix_z                       0.279 0.087      0.109    0.279      0.450
## ln_roads_km_c                  0.421 0.114      0.199    0.421      0.646
## roads_prop_highway_arterial_z  0.619 0.082      0.459    0.619      0.781
## ale_index_z                    0.488 0.178      0.140    0.488      0.837
## canbics_index_z                0.025 0.239     -0.445    0.026      0.494
## population_100_c               0.059 0.019      0.023    0.059      0.097
##                                 mode kld
## (Intercept)                   -0.208   0
## vandix_z                       0.279   0
## ln_roads_km_c                  0.421   0
## roads_prop_highway_arterial_z  0.619   0
## ale_index_z                    0.488   0
## canbics_index_z                0.026   0
## population_100_c               0.059   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.044 0.180      0.732    1.030      1.440 1.002
## Phi for ui       0.675 0.131      0.392    0.688      0.892 0.719
## 
## Deviance Information Criterion (DIC) ...............: 1271.74
## Deviance Information Criterion (DIC, saturated) ....: 620.13
## Effective number of parameters .....................: 173.99
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1264.84
## Effective number of parameters .................: 127.02
## 
## Marginal log-Likelihood:  -482.97 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
fraser_valley_results <- bind_rows(fraser_valley_all %>% mutate(Region = "Fraser Valley",Outcome = "All Injury Claims"),
                         fraser_valley_cyclist %>% mutate(Region = "Fraser Valley",Outcome = "Cyclist Injury Claims"),
                         fraser_valley_pedestrian %>% mutate(Region = "Fraser Valley",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

North Central

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

North Central

n_claims

359

31,422.0

87.5

180.5

0.0

1,906.0

0

22

6.1

n_casualty_claims

4,490.0

12.5

25.2

0.0

228.0

0

57

15.9

n_cyclist_claims

107.0

0.3

1.0

0.0

11.0

0

296

82.5

n_cyclist_casualty_claims

54.0

0.2

0.5

0.0

5.0

0

323

90.0

n_pedestrian_claims

241.0

0.7

2.0

0.0

27.0

0

249

69.4

n_pedestrian_casualty_claims

203.0

0.6

1.6

0.0

21.0

0

262

73.0

population

182,818.0

513.5

284.6

0.0

2,479.0

3

11

3.1

total_roads_km

9,989.6

27.8

53.9

0.0

398.4

0

9

2.5

roads_prop_highway_arterial

45.4

0.1

0.1

0.0

0.7

9

132

36.8

no_highschool_prevalance

74.8

0.2

0.1

0.1

0.8

30

0

0.0

unemployment_rate

3,459.5

10.5

6.8

0.0

60.0

30

14

3.9

hh_avg_income

24,402,848.0

77,716.1

21,870.6

34,760.0

201,050.0

45

0

0.0

participation_rate

22,327.9

67.9

9.3

35.7

89.8

30

0

0.0

university_degree_prevalance

149.0

0.5

0.1

0.0

0.7

30

1

0.3

lone_parent_fam_prevalence

60.1

0.2

0.1

0.0

0.6

29

0

0.0

home_owner_prevalence

233.8

0.7

0.2

0.0

1.0

30

4

1.1

vandix

156.6

0.5

0.7

-1.0

3.4

45

0

0.0

ale_index

-387.8

-1.3

0.6

-2.1

0.0

65

0

0.0

canbics_index

183.2

0.6

1.0

0.0

4.9

44

179

49.9

Maps

claims <- ggplot() + 
  geom_sf(data = north_central_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = north_central_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = north_central_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

north_central_sp <- as(north_central_sf, "Spatial")
north_central_sp$ui <- 1:nrow(north_central_sp@data)

coords <- coordinates(north_central_sp)
north_central_nb <- poly2nb(north_central_sp, queen = TRUE)

north_central_nb
## Neighbour list object:
## Number of regions: 359 
## Number of nonzero links: 1958 
## Percentage nonzero weights: 1.519231 
## Average number of links: 5.454039
#assign nearest neighbour for no links

north_central_nb <- assign_nearest_neighbors(north_central_nb,north_central_sp)

plot(north_central_sp, border = grey(0.5))
plot(north_central_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(north_central_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(north_central_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  north_central_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 7.6708, p-value = 8.548e-15
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.241279277      -0.002793296       0.001012416
cyc_cc_mi <- moran.test(north_central_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  north_central_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 3.7022, p-value = 0.0001069
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.116363461      -0.002793296       0.001035898
pd_cc_mi <- moran.test(north_central_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  north_central_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 6.5559, p-value = 2.766e-11
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1938375579     -0.0027932961      0.0008995886

BYM2 Models

nb2INLA("north_central.adj", north_central_nb)
g <- inla.read.graph(filename = "north_central.adj")


# Model Set 1: Total Casualty Claim Crashes

north_central_models_1 <- cma_models(north_central_sp@data, "North Central", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
north_central_all <- clean_fixed_effects(north_central_models_1$fixed_effects)

map(north_central_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.7, Running = 0.691, Post = 0.135, Total = 15.6 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 1.324 0.081      1.162    1.325      1.481 1.325   0
## vandix_z    0.191 0.076      0.042    0.191      0.340 0.191   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.296 0.044      0.218    0.293      0.391 0.287
## Phi for ui       0.788 0.062      0.650    0.794      0.890 0.807
## 
## Deviance Information Criterion (DIC) ...............: 1843.35
## Deviance Information Criterion (DIC, saturated) ....: 690.46
## Effective number of parameters .....................: 301.98
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1868.79
## Effective number of parameters .................: 223.04
## 
## Marginal log-Likelihood:  -1102.66 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.2, Running = 0.685, Post = 0.13, Total = 15 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   0.531 0.136      0.261    0.533      0.795 0.533   0
## vandix_z      0.272 0.072      0.132    0.272      0.413 0.272   0
## ln_roads_km_c 0.698 0.090      0.521    0.698      0.876 0.698   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.229 0.031      0.172    0.228      0.294 0.227
## Phi for ui       0.950 0.031      0.871    0.957      0.990 0.972
## 
## Deviance Information Criterion (DIC) ...............: 1801.00
## Deviance Information Criterion (DIC, saturated) ....: 648.12
## Effective number of parameters .....................: 283.39
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1797.13
## Effective number of parameters .................: 193.56
## 
## Marginal log-Likelihood:  -1079.96 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.1, Running = 0.755, Post = 0.156, Total = 15 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                    0.969 0.220      0.535    0.970      1.397
## vandix_z                       0.201 0.067      0.071    0.201      0.332
## ln_roads_km_c                  0.542 0.087      0.371    0.543      0.713
## roads_prop_highway_arterial_z  0.616 0.086      0.447    0.616      0.784
## ale_index_z                   -0.263 0.293     -0.837   -0.263      0.313
## canbics_index_z                0.384 0.311     -0.228    0.384      0.994
## population_100_c               0.034 0.021     -0.008    0.034      0.076
##                                 mode kld
## (Intercept)                    0.970   0
## vandix_z                       0.201   0
## ln_roads_km_c                  0.543   0
## roads_prop_highway_arterial_z  0.616   0
## ale_index_z                   -0.263   0
## canbics_index_z                0.384   0
## population_100_c               0.034   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.269 0.035      0.204    0.267      0.341 0.267
## Phi for ui       0.962 0.028      0.890    0.969      0.994 0.983
## 
## Deviance Information Criterion (DIC) ...............: 1795.38
## Deviance Information Criterion (DIC, saturated) ....: 642.50
## Effective number of parameters .....................: 274.27
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1785.78
## Effective number of parameters .................: 185.29
## 
## Marginal log-Likelihood:  -1076.84 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

north_central_models_2 <- cma_models(north_central_sp@data,"North Central","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
north_central_cyclist <- clean_fixed_effects(north_central_models_2$fixed_effects)

map(north_central_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.1, Running = 0.734, Post = 0.119, Total = 15 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -3.130 0.345     -3.906    -3.11     -2.501 -3.117   0
## vandix_z     0.311 0.138      0.043     0.31      0.584  0.310   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.649 0.201      0.347    0.618      1.131 0.559
## Phi for ui       0.158 0.109      0.024    0.131      0.438 0.071
## 
## Deviance Information Criterion (DIC) ...............: 277.05
## Deviance Information Criterion (DIC, saturated) ....: 195.66
## Effective number of parameters .....................: 53.64
## 
## Watanabe-Akaike information criterion (WAIC) ...: 281.36
## Effective number of parameters .................: 43.50
## 
## Marginal log-Likelihood:  -75.31 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14, Running = 0.743, Post = 0.124, Total = 14.9 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -3.335 0.401     -4.194   -3.317     -2.595 -3.320   0
## vandix_z       0.350 0.145      0.069    0.349      0.637  0.349   0
## ln_roads_km_c  0.160 0.162     -0.154    0.159      0.483  0.159   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.637 0.198      0.340    0.607      1.113 0.549
## Phi for ui       0.199 0.129      0.032    0.169      0.522 0.097
## 
## Deviance Information Criterion (DIC) ...............: 276.34
## Deviance Information Criterion (DIC, saturated) ....: 194.96
## Effective number of parameters .....................: 52.45
## 
## Watanabe-Akaike information criterion (WAIC) ...: 280.80
## Effective number of parameters .................: 43.03
## 
## Marginal log-Likelihood:  -80.11 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.2, Running = 0.783, Post = 0.136, Total = 15.1 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -2.521 0.500     -3.542   -2.507     -1.578
## vandix_z                       0.356 0.148      0.068    0.354      0.651
## ln_roads_km_c                  0.211 0.181     -0.141    0.210      0.570
## roads_prop_highway_arterial_z  0.504 0.199      0.115    0.503      0.898
## ale_index_z                   -0.474 0.713     -1.872   -0.475      0.929
## canbics_index_z                2.028 0.574      0.911    2.024      3.168
## population_100_c              -0.047 0.077     -0.197   -0.047      0.104
##                                 mode kld
## (Intercept)                   -2.508   0
## vandix_z                       0.354   0
## ln_roads_km_c                  0.210   0
## roads_prop_highway_arterial_z  0.503   0
## ale_index_z                   -0.475   0
## canbics_index_z                2.024   0
## population_100_c              -0.047   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.833 0.308      0.399    0.778      1.593 0.677
## Phi for ui       0.306 0.166      0.063    0.280      0.682 0.196
## 
## Deviance Information Criterion (DIC) ...............: 272.87
## Deviance Information Criterion (DIC, saturated) ....: 191.48
## Effective number of parameters .....................: 43.18
## 
## Watanabe-Akaike information criterion (WAIC) ...: 280.19
## Effective number of parameters .................: 39.40
## 
## Marginal log-Likelihood:  -89.42 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

north_central_models_3 <- cma_models(north_central_sp@data,"North Central","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
north_central_pedestrian <- clean_fixed_effects(north_central_models_3$fixed_effects)

map(north_central_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.1, Running = 0.72, Post = 0.122, Total = 15 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.929 0.196     -2.344   -1.921     -1.566 -1.922   0
## vandix_z     0.382 0.105      0.177    0.382      0.589  0.382   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.548 0.124      0.346    0.534      0.833 0.506
## Phi for ui       0.359 0.152      0.106    0.346      0.681 0.308
## 
## Deviance Information Criterion (DIC) ...............: 591.43
## Deviance Information Criterion (DIC, saturated) ....: 358.45
## Effective number of parameters .....................: 115.31
## 
## Watanabe-Akaike information criterion (WAIC) ...: 593.68
## Effective number of parameters .................: 84.99
## 
## Marginal log-Likelihood:  -253.06 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.1, Running = 0.752, Post = 0.137, Total = 15 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -2.058 0.246     -2.566   -2.050     -1.598 -2.051   0
## vandix_z       0.392 0.107      0.182    0.392      0.604  0.392   0
## ln_roads_km_c  0.120 0.141     -0.147    0.116      0.409  0.116   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.512 0.132      0.302    0.495      0.818 0.463
## Phi for ui       0.451 0.180      0.131    0.445      0.801 0.425
## 
## Deviance Information Criterion (DIC) ...............: 588.55
## Deviance Information Criterion (DIC, saturated) ....: 355.56
## Effective number of parameters .....................: 112.61
## 
## Watanabe-Akaike information criterion (WAIC) ...: 590.64
## Effective number of parameters .................: 83.39
## 
## Marginal log-Likelihood:  -258.22 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.1, Running = 0.79, Post = 0.154, Total = 15.1 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -2.107 0.444     -3.013   -2.095     -1.271
## vandix_z                       0.335 0.107      0.126    0.335      0.548
## ln_roads_km_c                 -0.063 0.162     -0.372   -0.066      0.265
## roads_prop_highway_arterial_z  0.650 0.146      0.364    0.650      0.939
## ale_index_z                   -1.238 0.631     -2.498   -1.230     -0.021
## canbics_index_z                0.642 0.503     -0.340    0.640      1.637
## population_100_c               0.022 0.044     -0.064    0.022      0.109
##                                 mode kld
## (Intercept)                   -2.096   0
## vandix_z                       0.335   0
## ln_roads_km_c                 -0.066   0
## roads_prop_highway_arterial_z  0.650   0
## ale_index_z                   -1.230   0
## canbics_index_z                0.640   0
## population_100_c               0.022   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.522 0.147      0.291    0.503      0.864 0.466
## Phi for ui       0.560 0.181      0.204    0.568      0.875 0.598
## 
## Deviance Information Criterion (DIC) ...............: 575.77
## Deviance Information Criterion (DIC, saturated) ....: 342.79
## Effective number of parameters .....................: 103.09
## 
## Watanabe-Akaike information criterion (WAIC) ...: 578.23
## Effective number of parameters .................: 77.80
## 
## Marginal log-Likelihood:  -266.22 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
north_central_results <- bind_rows(north_central_all %>% mutate(Region = "North Central",Outcome = "All Injury Claims"),
                         north_central_cyclist %>% mutate(Region = "North Central",Outcome = "Cyclist Injury Claims"),
                         north_central_pedestrian %>% mutate(Region = "North Central",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

Kamloops-Salmon Arm

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Kamloops-Salmon Arm

n_claims

204

24,765.0

121.4

244.3

0.0

2,286.0

0

5

2.5

n_casualty_claims

4,069.0

19.9

39.5

0.0

293.0

0

19

9.3

n_cyclist_claims

143.0

0.7

1.7

0.0

14.0

0

138

67.6

n_cyclist_casualty_claims

89.0

0.4

1.1

0.0

10.0

0

152

74.5

n_pedestrian_claims

236.0

1.2

2.9

0.0

21.0

0

130

63.7

n_pedestrian_casualty_claims

183.0

0.9

2.3

0.0

17.0

0

137

67.2

population

133,847.0

656.1

348.0

0.0

2,310.0

0

2

1.0

total_roads_km

2,744.5

13.5

25.3

0.5

217.0

0

0

0.0

roads_prop_highway_arterial

49.9

0.2

0.2

0.0

0.7

0

32

15.7

no_highschool_prevalance

29.4

0.2

0.1

0.0

0.4

31

0

0.0

unemployment_rate

1,439.0

8.3

5.3

0.0

37.5

31

12

5.9

hh_avg_income

12,008,631.0

71,479.9

20,375.9

29,162.0

137,275.0

36

0

0.0

participation_rate

10,775.2

62.3

9.9

30.6

90.0

31

0

0.0

university_degree_prevalance

88.4

0.5

0.1

0.3

0.7

31

0

0.0

lone_parent_fam_prevalence

31.0

0.2

0.1

0.0

0.5

30

0

0.0

home_owner_prevalence

128.1

0.7

0.2

0.1

1.0

31

0

0.0

vandix

30.5

0.2

0.7

-1.3

2.4

36

0

0.0

ale_index

-116.3

-0.7

1.2

-2.1

3.3

41

1

0.5

canbics_index

227.8

1.4

1.5

0.0

5.6

41

40

19.6

Maps

claims <- ggplot() + 
  geom_sf(data = interior_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = interior_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = interior_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

interior_sp <- as(interior_sf, "Spatial")
interior_sp$ui <- 1:nrow(interior_sp@data)

coords <- coordinates(interior_sp)
interior_nb <- poly2nb(interior_sp, queen = TRUE)
## Warning in poly2nb(interior_sp, queen = TRUE): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
interior_nb
## Neighbour list object:
## Number of regions: 204 
## Number of nonzero links: 1098 
## Percentage nonzero weights: 2.638408 
## Average number of links: 5.382353 
## 3 disjoint connected subgraphs
#assign nearest neighbour for no links

interior_nb <- assign_nearest_neighbors(interior_nb,interior_sp)

plot(interior_sp, border = grey(0.5))
plot(interior_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(interior_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(interior_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  interior_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 6.3114, p-value = 1.383e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.261379586      -0.004926108       0.001780360
cyc_cc_mi <- moran.test(interior_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  interior_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 2.4815, p-value = 0.006542
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.096406801      -0.004926108       0.001667545
pd_cc_mi <- moran.test(interior_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  interior_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 4.6131, p-value = 1.983e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.190275369      -0.004926108       0.001790493

BYM2 Models

nb2INLA("interior.adj", interior_nb)
g <- inla.read.graph(filename = "interior.adj")


# Model Set 1: Total Casualty Claim Crashes

interior_models_1 <- cma_models(interior_sp@data, "Kamloops-Salmon Arm", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
interior_all <- clean_fixed_effects(interior_models_1$fixed_effects)

map(interior_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.6, Running = 0.538, Post = 0.114, Total = 15.3 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 1.953 0.082      1.790    1.954      2.112 1.954   0
## vandix_z    0.250 0.115      0.024    0.251      0.476 0.251   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.647 0.083      0.498    0.642      0.825 0.632
## Phi for ui       0.421 0.107      0.223    0.417      0.639 0.409
## 
## Deviance Information Criterion (DIC) ...............: 1176.14
## Deviance Information Criterion (DIC, saturated) ....: 411.35
## Effective number of parameters .....................: 182.29
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1181.66
## Effective number of parameters .................: 129.25
## 
## Marginal log-Likelihood:  -647.01 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.6, Running = 0.519, Post = 0.127, Total = 15.3 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   1.431 0.099      1.234    1.431      1.624 1.431   0
## vandix_z      0.258 0.104      0.055    0.258      0.461 0.258   0
## ln_roads_km_c 0.850 0.106      0.644    0.849      1.058 0.849   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.776 0.104      0.590     0.77      0.999 0.759
## Phi for ui       0.625 0.097      0.425     0.63      0.802 0.637
## 
## Deviance Information Criterion (DIC) ...............: 1158.35
## Deviance Information Criterion (DIC, saturated) ....: 393.55
## Effective number of parameters .....................: 172.37
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1151.76
## Effective number of parameters .................: 115.90
## 
## Marginal log-Likelihood:  -622.71 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.5, Running = 0.553, Post = 0.145, Total = 15.2 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                    1.132 0.119      0.898    1.133      1.366
## vandix_z                       0.208 0.085      0.042    0.208      0.375
## ln_roads_km_c                  0.491 0.096      0.302    0.491      0.681
## roads_prop_highway_arterial_z  0.534 0.056      0.424    0.534      0.645
## ale_index_z                    0.091 0.257     -0.413    0.091      0.594
## canbics_index_z               -0.712 0.249     -1.202   -0.712     -0.224
## population_100_c               0.099 0.019      0.062    0.099      0.136
##                                 mode kld
## (Intercept)                    1.133   0
## vandix_z                       0.208   0
## ln_roads_km_c                  0.491   0
## roads_prop_highway_arterial_z  0.534   0
## ale_index_z                    0.091   0
## canbics_index_z               -0.712   0
## population_100_c               0.099   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.176 0.170      0.871    1.166      1.540 1.149
## Phi for ui       0.819 0.085      0.620    0.832      0.948 0.862
## 
## Deviance Information Criterion (DIC) ...............: 1138.16
## Deviance Information Criterion (DIC, saturated) ....: 373.37
## Effective number of parameters .....................: 154.77
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1123.45
## Effective number of parameters .................: 100.52
## 
## Marginal log-Likelihood:  -598.04 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

interior_models_2 <- cma_models(interior_sp@data,"Kamloops-Salmon Arm","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
interior_cyclist <- clean_fixed_effects(interior_models_2$fixed_effects)

map(interior_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.5, Running = 0.524, Post = 0.131, Total = 15.2 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.615 0.208     -2.054   -1.608     -1.228 -1.609   0
## vandix_z     0.372 0.154      0.064    0.373      0.668  0.374   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.061 0.334      0.563    1.009      1.863 0.911
## Phi for ui       0.209 0.131      0.034    0.182      0.528 0.106
## 
## Deviance Information Criterion (DIC) ...............: 323.23
## Deviance Information Criterion (DIC, saturated) ....: 202.50
## Effective number of parameters .....................: 50.21
## 
## Watanabe-Akaike information criterion (WAIC) ...: 325.30
## Effective number of parameters .................: 40.32
## 
## Marginal log-Likelihood:  -58.43 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.5, Running = 0.53, Post = 0.147, Total = 15.1 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -2.040 0.278     -2.618   -2.030     -1.523 -2.031   0
## vandix_z       0.394 0.163      0.069    0.395      0.709  0.395   0
## ln_roads_km_c  0.586 0.201      0.200    0.583      0.990  0.583   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.066 0.349      0.551     1.01      1.909 0.906
## Phi for ui       0.519 0.193      0.160     0.52      0.871 0.514
## 
## Deviance Information Criterion (DIC) ...............: 316.59
## Deviance Information Criterion (DIC, saturated) ....: 195.85
## Effective number of parameters .....................: 45.63
## 
## Watanabe-Akaike information criterion (WAIC) ...: 318.57
## Effective number of parameters .................: 37.17
## 
## Marginal log-Likelihood:  -59.22 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.5, Running = 0.568, Post = 0.122, Total = 15.2 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -1.919 0.324     -2.590   -1.907     -1.318
## vandix_z                       0.391 0.177      0.034    0.395      0.726
## ln_roads_km_c                  0.195 0.216     -0.224    0.193      0.624
## roads_prop_highway_arterial_z  0.395 0.138      0.125    0.394      0.668
## ale_index_z                    0.630 0.488     -0.340    0.635      1.576
## canbics_index_z               -0.284 0.489     -1.239   -0.286      0.683
## population_100_c               0.141 0.038      0.067    0.141      0.217
##                                 mode kld
## (Intercept)                   -1.907   0
## vandix_z                       0.395   0
## ln_roads_km_c                  0.193   0
## roads_prop_highway_arterial_z  0.394   0
## ale_index_z                    0.635   0
## canbics_index_z               -0.286   0
## population_100_c               0.141   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.847 0.897      0.722    1.648       4.15 1.324
## Phi for ui       0.619 0.218      0.170    0.646       0.95 0.803
## 
## Deviance Information Criterion (DIC) ...............: 308.25
## Deviance Information Criterion (DIC, saturated) ....: 187.51
## Effective number of parameters .....................: 34.28
## 
## Watanabe-Akaike information criterion (WAIC) ...: 311.20
## Effective number of parameters .................: 30.41
## 
## Marginal log-Likelihood:  -66.20 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

interior_models_3 <- cma_models(interior_sp@data,"Kamloops-Salmon Arm","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
interior_pedestrian <- clean_fixed_effects(interior_models_3$fixed_effects)

map(interior_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.5, Running = 0.532, Post = 0.116, Total = 15.1 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.463 0.201     -1.886   -1.455     -1.091 -1.456   0
## vandix_z     0.654 0.159      0.343    0.653      0.970  0.653   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.626 0.133      0.404    0.612      0.925 0.586
## Phi for ui       0.267 0.124      0.075    0.250      0.549 0.207
## 
## Deviance Information Criterion (DIC) ...............: 415.59
## Deviance Information Criterion (DIC, saturated) ....: 244.14
## Effective number of parameters .....................: 81.09
## 
## Watanabe-Akaike information criterion (WAIC) ...: 423.72
## Effective number of parameters .................: 63.39
## 
## Marginal log-Likelihood:  -121.84 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.6, Running = 0.53, Post = 0.149, Total = 15.2 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -1.851 0.243     -2.350   -1.843     -1.394 -1.844   0
## vandix_z       0.693 0.162      0.377    0.692      1.012  0.692   0
## ln_roads_km_c  0.620 0.181      0.269    0.619      0.980  0.619   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.676 0.152      0.426     0.66      1.020 0.628
## Phi for ui       0.474 0.155      0.190     0.47      0.778 0.454
## 
## Deviance Information Criterion (DIC) ...............: 408.87
## Deviance Information Criterion (DIC, saturated) ....: 237.41
## Effective number of parameters .....................: 73.21
## 
## Watanabe-Akaike information criterion (WAIC) ...: 414.87
## Effective number of parameters .................: 57.58
## 
## Marginal log-Likelihood:  -121.42 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.5, Running = 0.564, Post = 0.123, Total = 15.2 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -1.579 0.290     -2.173   -1.570     -1.035
## vandix_z                       0.665 0.153      0.362    0.665      0.966
## ln_roads_km_c                  0.248 0.189     -0.120    0.247      0.622
## roads_prop_highway_arterial_z  0.389 0.123      0.148    0.388      0.633
## ale_index_z                    1.566 0.450      0.675    1.568      2.446
## canbics_index_z               -0.983 0.452     -1.881   -0.980     -0.103
## population_100_c               0.111 0.037      0.039    0.111      0.184
##                                 mode kld
## (Intercept)                   -1.570   0
## vandix_z                       0.665   0
## ln_roads_km_c                  0.247   0
## roads_prop_highway_arterial_z  0.388   0
## ale_index_z                    1.568   0
## canbics_index_z               -0.980   0
## population_100_c               0.111   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.083 0.313      0.602    1.038      1.823 0.953
## Phi for ui       0.373 0.197      0.064    0.352      0.783 0.253
## 
## Deviance Information Criterion (DIC) ...............: 400.81
## Deviance Information Criterion (DIC, saturated) ....: 229.36
## Effective number of parameters .....................: 61.39
## 
## Watanabe-Akaike information criterion (WAIC) ...: 404.56
## Effective number of parameters .................: 49.31
## 
## Marginal log-Likelihood:  -125.78 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
interior_results <- bind_rows(interior_all %>% mutate(Region = "Kamloops-Salmon Arm",Outcome = "All Injury Claims"),
                         interior_cyclist %>% mutate(Region = "Kamloops-Salmon Arm",Outcome = "Cyclist Injury Claims"),
                         interior_pedestrian %>% mutate(Region = "Kamloops-Salmon Arm",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

Southeast

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Southeast

n_claims

110

9,978.0

90.7

170.1

0.0

1,235.0

0

5

4.5

n_casualty_claims

1,037.0

9.4

17.8

0.0

121.0

0

14

12.7

n_cyclist_claims

40.0

0.4

0.8

0.0

4.0

0

87

79.1

n_cyclist_casualty_claims

24.0

0.2

0.6

0.0

4.0

0

94

85.5

n_pedestrian_claims

62.0

0.6

1.3

0.0

8.0

0

79

71.8

n_pedestrian_casualty_claims

51.0

0.5

1.2

0.0

8.0

0

85

77.3

population

60,427.0

549.3

214.3

0.0

1,482.0

0

4

3.6

total_roads_km

1,598.6

14.5

23.7

0.0

202.4

0

2

1.8

roads_prop_highway_arterial

13.9

0.1

0.2

0.0

0.5

2

49

44.5

no_highschool_prevalance

16.2

0.2

0.1

0.0

0.3

6

0

0.0

unemployment_rate

838.5

8.1

4.4

0.0

17.5

6

9

8.2

hh_avg_income

7,149,722.0

68,747.3

16,755.0

34,159.0

122,980.0

6

0

0.0

participation_rate

6,473.1

62.2

8.5

40.9

80.0

6

0

0.0

university_degree_prevalance

57.1

0.5

0.1

0.3

0.7

6

0

0.0

lone_parent_fam_prevalence

17.1

0.2

0.1

0.0

0.4

6

0

0.0

home_owner_prevalence

77.2

0.7

0.2

0.0

1.0

6

1

0.9

vandix

8.2

0.1

0.6

-0.9

1.6

6

0

0.0

ale_index

-125.7

-1.2

0.7

-2.1

0.4

9

0

0.0

canbics_index

32.4

0.3

0.5

0.0

2.0

9

64

58.2

Maps

claims <- ggplot() + 
  geom_sf(data = southeast_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = southeast_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = southeast_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

southeast_sp <- as(southeast_sf, "Spatial")
southeast_sp$ui <- 1:nrow(southeast_sp@data)

coords <- coordinates(southeast_sp)
southeast_nb <- poly2nb(southeast_sp, queen = TRUE)
## Warning in poly2nb(southeast_sp, queen = TRUE): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
southeast_nb
## Neighbour list object:
## Number of regions: 110 
## Number of nonzero links: 552 
## Percentage nonzero weights: 4.561983 
## Average number of links: 5.018182 
## 3 disjoint connected subgraphs
#assign nearest neighbour for no links

southeast_nb <- assign_nearest_neighbors(southeast_nb,southeast_sp)

plot(southeast_sp, border = grey(0.5))
plot(southeast_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(southeast_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(southeast_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  southeast_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 2.7069, p-value = 0.003396
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.146644926      -0.009174312       0.003313622
cyc_cc_mi <- moran.test(southeast_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  southeast_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 0.65675, p-value = 0.2557
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.029151043      -0.009174312       0.003405459
pd_cc_mi <- moran.test(southeast_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  southeast_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 2.1446, p-value = 0.01599
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.114119689      -0.009174312       0.003305202

BYM2 Models

nb2INLA("southeast.adj", southeast_nb)
g <- inla.read.graph(filename = "southeast.adj")


# Model Set 1: Total Casualty Claim Crashes

southeast_models_1 <- cma_models(southeast_sp@data, "Southeast", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
southeast_all <- clean_fixed_effects(southeast_models_1$fixed_effects)

map(southeast_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.4, Running = 0.42, Post = 0.14, Total = 14.9 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 1.270 0.111      1.046    1.271      1.484 1.271   0
## vandix_z    0.352 0.150      0.057    0.352      0.648 0.352   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.750 0.128      0.528    0.740      1.029 0.721
## Phi for ui       0.396 0.141      0.151    0.387      0.689 0.361
## 
## Deviance Information Criterion (DIC) ...............: 543.48
## Deviance Information Criterion (DIC, saturated) ....: 207.17
## Effective number of parameters .....................: 89.95
## 
## Watanabe-Akaike information criterion (WAIC) ...: 547.74
## Effective number of parameters .................: 64.68
## 
## Marginal log-Likelihood:  -266.28 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.4, Running = 0.434, Post = 0.15, Total = 14.9 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   0.651 0.166      0.321    0.652      0.973 0.652   0
## vandix_z      0.386 0.141      0.111    0.385      0.664 0.385   0
## ln_roads_km_c 0.779 0.157      0.474    0.778      1.092 0.778   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.827 0.147      0.571    0.815      1.150 0.795
## Phi for ui       0.691 0.145      0.377    0.707      0.923 0.754
## 
## Deviance Information Criterion (DIC) ...............: 530.90
## Deviance Information Criterion (DIC, saturated) ....: 194.60
## Effective number of parameters .....................: 83.84
## 
## Watanabe-Akaike information criterion (WAIC) ...: 528.31
## Effective number of parameters .................: 56.64
## 
## Marginal log-Likelihood:  -259.52 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.5, Running = 0.524, Post = 0.126, Total = 15.1 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                    1.564 0.397      0.776    1.566      2.338
## vandix_z                       0.350 0.131      0.094    0.349      0.609
## ln_roads_km_c                  0.490 0.179      0.142    0.489      0.844
## roads_prop_highway_arterial_z  0.627 0.150      0.331    0.627      0.921
## ale_index_z                    1.634 0.820      0.032    1.631      3.253
## canbics_index_z               -0.712 0.851     -2.390   -0.711      0.955
## population_100_c               0.152 0.057      0.040    0.151      0.265
##                                 mode kld
## (Intercept)                    1.566   0
## vandix_z                       0.349   0
## ln_roads_km_c                  0.489   0
## roads_prop_highway_arterial_z  0.627   0
## ale_index_z                    1.631   0
## canbics_index_z               -0.711   0
## population_100_c               0.151   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.994 0.183      0.677    0.979      1.397 0.952
## Phi for ui       0.771 0.134      0.457    0.795      0.961 0.870
## 
## Deviance Information Criterion (DIC) ...............: 524.87
## Deviance Information Criterion (DIC, saturated) ....: 188.56
## Effective number of parameters .....................: 79.56
## 
## Watanabe-Akaike information criterion (WAIC) ...: 518.83
## Effective number of parameters .................: 52.02
## 
## Marginal log-Likelihood:  -267.01 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

southeast_models_2 <- cma_models(southeast_sp@data,"Southeast","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
southeast_cyclist <- clean_fixed_effects(southeast_models_2$fixed_effects)

map(southeast_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.4, Running = 0.4, Post = 0.147, Total = 15 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.997 0.353     -2.783   -1.966     -1.387 -1.913   0
## vandix_z     0.597 0.230      0.156    0.593      1.061  0.593   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                    mean     sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 11.145 30.615      0.759    4.049     67.116 1.404
## Phi for ui        0.277  0.217      0.018    0.218      0.796 0.048
## 
## Deviance Information Criterion (DIC) ...............: 125.20
## Deviance Information Criterion (DIC, saturated) ....: 88.73
## Effective number of parameters .....................: 9.52
## 
## Watanabe-Akaike information criterion (WAIC) ...: 133.43
## Effective number of parameters .................: 15.36
## 
## Marginal log-Likelihood:  11.87 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.4, Running = 0.424, Post = 0.126, Total = 14.9 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -2.450 0.436     -3.362   -2.433     -1.643 -2.434   0
## vandix_z       0.755 0.249      0.271    0.753      1.249  0.753   0
## ln_roads_km_c  0.468 0.243     -0.008    0.467      0.948  0.467   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                    mean      sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 60.416 409.765      0.709    7.817     408.21 1.435
## Phi for ui        0.357   0.246      0.027    0.308       0.88 0.078
## 
## Deviance Information Criterion (DIC) ...............: 123.94
## Deviance Information Criterion (DIC, saturated) ....: 87.47
## Effective number of parameters .....................: 9.67
## 
## Watanabe-Akaike information criterion (WAIC) ...: 128.56
## Effective number of parameters .................: 12.56
## 
## Marginal log-Likelihood:  9.28 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.5, Running = 0.425, Post = 0.128, Total = 15.1 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -0.882 0.610     -2.094   -0.877      0.301
## vandix_z                       0.828 0.252      0.337    0.828      1.325
## ln_roads_km_c                  0.202 0.376     -0.535    0.202      0.939
## roads_prop_highway_arterial_z  0.331 0.334     -0.323    0.330      0.989
## ale_index_z                    1.940 1.360     -0.715    1.936      4.622
## canbics_index_z                0.074 1.281     -2.441    0.075      2.586
## population_100_c               0.288 0.094      0.105    0.288      0.473
##                                 mode kld
## (Intercept)                   -0.878   0
## vandix_z                       0.828   0
## ln_roads_km_c                  0.202   0
## roads_prop_highway_arterial_z  0.330   0
## ale_index_z                    1.937   0
## canbics_index_z                0.075   0
## population_100_c               0.288   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                      mean       sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 5986.857 1.50e+05      0.933   66.570   2.28e+04 0.970
## Phi for ui          0.353 2.63e-01      0.018    0.293   9.06e-01 0.041
## 
## Deviance Information Criterion (DIC) ...............: 120.86
## Deviance Information Criterion (DIC, saturated) ....: 84.39
## Effective number of parameters .....................: 7.95
## 
## Watanabe-Akaike information criterion (WAIC) ...: 124.65
## Effective number of parameters .................: 10.26
## 
## Marginal log-Likelihood:  -0.811 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

southeast_models_3 <- cma_models(southeast_sp@data,"Southeast","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
southeast_pedestrian <- clean_fixed_effects(southeast_models_3$fixed_effects)

map(southeast_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.4, Running = 0.443, Post = 0.112, Total = 14.9 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.868 0.318     -2.559   -1.851     -1.289 -1.807   0
## vandix_z     0.545 0.240      0.075    0.543      1.021  0.543   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.728 0.242      0.368    0.689      1.310 0.618
## Phi for ui       0.165 0.143      0.011    0.122      0.546 0.031
## 
## Deviance Information Criterion (DIC) ...............: 167.70
## Deviance Information Criterion (DIC, saturated) ....: 107.52
## Effective number of parameters .....................: 33.05
## 
## Watanabe-Akaike information criterion (WAIC) ...: 171.01
## Effective number of parameters .................: 26.27
## 
## Marginal log-Likelihood:  -23.60 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14.5, Running = 0.408, Post = 0.157, Total = 15 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -2.101 0.378     -2.898   -2.085     -1.403 -2.087   0
## vandix_z       0.597 0.247      0.113    0.597      1.085  0.597   0
## ln_roads_km_c  0.295 0.257     -0.211    0.295      0.801  0.295   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.753 0.260      0.372    0.711      1.382 0.633
## Phi for ui       0.204 0.167      0.014    0.156      0.634 0.039
## 
## Deviance Information Criterion (DIC) ...............: 168.93
## Deviance Information Criterion (DIC, saturated) ....: 108.74
## Effective number of parameters .....................: 32.75
## 
## Watanabe-Akaike information criterion (WAIC) ...: 173.11
## Effective number of parameters .................: 26.76
## 
## Marginal log-Likelihood:  -27.78 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14.6, Running = 0.475, Post = 0.147, Total = 15.2 
## Fixed effects:
##                                 mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -1.384 0.793     -2.998   -1.366      0.125
## vandix_z                       0.606 0.253      0.111    0.605      1.106
## ln_roads_km_c                 -0.070 0.387     -0.840   -0.067      0.685
## roads_prop_highway_arterial_z  0.307 0.332     -0.349    0.308      0.958
## ale_index_z                    0.615 1.423     -2.208    0.622      3.394
## canbics_index_z               -0.167 1.444     -2.983   -0.177      2.704
## population_100_c               0.233 0.115      0.010    0.232      0.464
##                                 mode kld
## (Intercept)                   -1.367   0
## vandix_z                       0.606   0
## ln_roads_km_c                 -0.067   0
## roads_prop_highway_arterial_z  0.308   0
## ale_index_z                    0.622   0
## canbics_index_z               -0.177   0
## population_100_c               0.232   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.736 0.269      0.349    0.690      1.392 0.606
## Phi for ui       0.227 0.180      0.016    0.177      0.679 0.045
## 
## Deviance Information Criterion (DIC) ...............: 170.16
## Deviance Information Criterion (DIC, saturated) ....: 109.97
## Effective number of parameters .....................: 33.93
## 
## Watanabe-Akaike information criterion (WAIC) ...: 177.67
## Effective number of parameters .................: 29.53
## 
## Marginal log-Likelihood:  -42.54 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
southeast_results <- bind_rows(southeast_all %>% mutate(Region = "Southeast",Outcome = "All Injury Claims"),
                         southeast_cyclist %>% mutate(Region = "Southeast",Outcome = "Cyclist Injury Claims"),
                         southeast_pedestrian %>% mutate(Region = "Southeast",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

Northwest

Descriptive Statistics

region

variable

n

sum

mean

sd

min

max

missing

n_zero

p_zero

Northwest

n_claims

78

4,596.0

58.9

105.9

0.0

746.0

0

5

6.4

n_casualty_claims

498.0

6.4

11.8

0.0

67.0

0

17

21.8

n_cyclist_claims

18.0

0.2

0.6

0.0

3.0

0

67

85.9

n_cyclist_casualty_claims

9.0

0.1

0.4

0.0

2.0

0

70

89.7

n_pedestrian_claims

51.0

0.7

1.2

0.0

5.0

0

54

69.2

n_pedestrian_casualty_claims

41.0

0.5

1.1

0.0

5.0

0

57

73.1

population

33,048.0

434.8

139.6

0.0

747.0

2

1

1.3

total_roads_km

1,124.4

14.4

29.7

0.0

151.7

0

3

3.8

roads_prop_highway_arterial

11.7

0.2

0.2

0.0

1.0

3

36

46.2

no_highschool_prevalance

18.6

0.2

0.1

0.1

0.6

3

0

0.0

unemployment_rate

934.5

12.5

9.3

0.0

50.0

3

4

5.1

hh_avg_income

5,356,607.0

74,397.3

16,951.4

39,451.0

115,632.0

6

0

0.0

participation_rate

5,081.6

67.8

7.7

36.4

81.4

3

0

0.0

university_degree_prevalance

34.3

0.5

0.1

0.2

0.6

3

0

0.0

lone_parent_fam_prevalence

14.7

0.2

0.1

0.0

0.4

3

0

0.0

home_owner_prevalence

52.5

0.7

0.2

0.2

1.0

3

0

0.0

vandix

48.8

0.7

0.9

-0.7

3.5

6

0

0.0

ale_index

-108.7

-1.5

0.5

-2.1

-0.5

5

0

0.0

canbics_index

0.0

0.0

0.0

0.0

0.0

4

74

94.9

Maps

claims <- ggplot() + 
  geom_sf(data = northwest_sf, aes(fill = n_casualty_claims,colour=n_casualty_claims)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Insurance Claims",
                     type = "aggregation", palette = "Earth", direction = -1) + 
  theme_void() + 
  ggtitle(
    "Vancouver"
  )


vandix <- ggplot() + 
  geom_sf(data = northwest_sf, aes(fill = vandix_z_c,colour=vandix_z_c)) +
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  scale_colour_carto_d(name = "VanDIX Score ",
                     type = "diverging", palette = "Earth", direction = -1) + 
  theme_void()


total_roads <- ggplot() + 
  geom_sf(data = northwest_sf, aes(fill = total_roads_km,colour=total_roads_km)) + 
  coord_sf(crs = "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs") + 
  scale_fill_carto_c(name = "Kilometres of Road",
                     type = "quantitative", palette = "Earth", direction = -1) + 
  scale_colour_carto_c(name = "Kilometres of Road",
                       type = "quantitative", palette = "Earth", direction = -1) + 
  theme_void()

cowplot::plot_grid(claims,vandix,total_roads,ncol=1)

Setting up Spatial Weights Matrix

#### Define Spatial Neighrbourhoods

northwest_sp <- as(northwest_sf, "Spatial")
northwest_sp$ui <- 1:nrow(northwest_sp@data)

coords <- coordinates(northwest_sp)
northwest_nb <- poly2nb(northwest_sp, queen = TRUE)

northwest_nb
## Neighbour list object:
## Number of regions: 78 
## Number of nonzero links: 406 
## Percentage nonzero weights: 6.673241 
## Average number of links: 5.205128
#assign nearest neighbour for no links

northwest_nb <- assign_nearest_neighbors(northwest_nb,northwest_sp)

plot(northwest_sp, border = grey(0.5))
plot(northwest_nb,
     coords = coords,
     add = TRUE, pch = 16, lwd = 2)

Spatial Dependency in each crash type

listw <- nb2listw(northwest_nb,zero.policy = TRUE)
all_cc_mi <- moran.test(northwest_sf$n_casualty_claims, listw,zero.policy = TRUE)
all_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  northwest_sf$n_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 4.5914, p-value = 2.201e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.290791318      -0.012987013       0.004377475
cyc_cc_mi <- moran.test(northwest_sf$n_cyclist_casualty_claims, listw,zero.policy = TRUE)
cyc_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  northwest_sf$n_cyclist_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 1.1465, p-value = 0.1258
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.064635043      -0.012987013       0.004583381
pd_cc_mi <- moran.test(northwest_sf$n_pedestrian_casualty_claims, listw,zero.policy = TRUE)
pd_cc_mi
## 
##  Moran I test under randomisation
## 
## data:  northwest_sf$n_pedestrian_casualty_claims  
## weights: listw    
## 
## Moran I statistic standard deviate = 4.4101, p-value = 5.165e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.298300116      -0.012987013       0.004982139

BYM2 Models

nb2INLA("northwest.adj", northwest_nb)
g <- inla.read.graph(filename = "northwest.adj")


# Model Set 1: Total Casualty Claim Crashes

northwest_models_1 <- cma_models(northwest_sp@data, "Northwest", "n_casualty_claims", "vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
northwest_all <- clean_fixed_effects(northwest_models_1$fixed_effects)

map(northwest_models_1$models, ~ summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.1, Running = 0.384, Post = 0.152, Total = 14.6 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.682 0.196      0.281    0.687      1.053 0.687   0
## vandix_z    0.281 0.122      0.042    0.280      0.521 0.280   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.632 0.158      0.381    0.612      0.999 0.574
## Phi for ui       0.367 0.202      0.054    0.345      0.785 0.226
## 
## Deviance Information Criterion (DIC) ...............: 354.83
## Deviance Information Criterion (DIC, saturated) ....: 151.31
## Effective number of parameters .....................: 61.55
## 
## Watanabe-Akaike information criterion (WAIC) ...: 364.65
## Effective number of parameters .................: 48.83
## 
## Marginal log-Likelihood:  -189.93 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14, Running = 0.388, Post = 0.138, Total = 14.5 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)   0.162 0.216     -0.275    0.166      0.573 0.166   0
## vandix_z      0.362 0.115      0.138    0.361      0.591 0.361   0
## ln_roads_km_c 0.780 0.160      0.469    0.780      1.096 0.780   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 0.559 0.144      0.323    0.543      0.885 0.515
## Phi for ui       0.810 0.148      0.430    0.851      0.983 0.950
## 
## Deviance Information Criterion (DIC) ...............: 338.75
## Deviance Information Criterion (DIC, saturated) ....: 135.24
## Effective number of parameters .....................: 54.97
## 
## Watanabe-Akaike information criterion (WAIC) ...: 340.10
## Effective number of parameters .................: 39.26
## 
## Marginal log-Likelihood:  -184.44 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14, Running = 0.418, Post = 0.123, Total = 14.5 
## Fixed effects:
##                                  mean    sd 0.025quant 0.5quant 0.975quant
## (Intercept)                    -5.270 8.504    -21.945   -5.270     11.410
## vandix_z                        0.222 0.107      0.013    0.222      0.433
## ln_roads_km_c                   0.423 0.284     -0.130    0.422      0.986
## roads_prop_highway_arterial_z   0.228 0.221     -0.209    0.228      0.663
## ale_index_z                     3.505 0.914      1.692    3.510      5.292
## canbics_index_z               -10.849 9.886    -30.222  -10.854      8.552
## population_100_c                0.372 0.123      0.130    0.373      0.614
##                                  mode kld
## (Intercept)                    -5.270   0
## vandix_z                        0.222   0
## ln_roads_km_c                   0.422   0
## roads_prop_highway_arterial_z   0.228   0
## ale_index_z                     3.510   0
## canbics_index_z               -10.854   0
## population_100_c                0.373   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.159 0.296      0.688    1.121      1.846 1.049
## Phi for ui       0.192 0.214      0.002    0.104      0.776 0.002
## 
## Deviance Information Criterion (DIC) ...............: 353.41
## Deviance Information Criterion (DIC, saturated) ....: 149.90
## Effective number of parameters .....................: 61.24
## 
## Watanabe-Akaike information criterion (WAIC) ...: 366.70
## Effective number of parameters .................: 50.50
## 
## Marginal log-Likelihood:  -186.55 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 2: Total Casualty Cyclist Claim Crashes

northwest_models_2 <- cma_models(northwest_sp@data,"Northwest","n_cyclist_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
northwest_cyclist <- clean_fixed_effects(northwest_models_2$fixed_effects)

map(northwest_models_2$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 13.9, Running = 0.419, Post = 0.133, Total = 14.5 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -2.583 0.466     -3.500   -2.582     -1.670 -2.582   0
## vandix_z     0.237 0.210     -0.173    0.237      0.649  0.237   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                      mean       sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1271.475 14399.03      2.365   86.105   7984.694 3.618
## Phi for ui          0.343     0.28      0.009    0.269      0.927 0.012
## 
## Deviance Information Criterion (DIC) ...............: 60.72
## Deviance Information Criterion (DIC, saturated) ....: 44.01
## Effective number of parameters .....................: 2.28
## 
## Watanabe-Akaike information criterion (WAIC) ...: 60.86
## Effective number of parameters .................: 2.37
## 
## Marginal log-Likelihood:  6.05 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 14, Running = 0.368, Post = 0.142, Total = 14.5 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -2.630 0.519     -3.651   -2.630     -1.615 -2.630   0
## vandix_z       0.246 0.217     -0.180    0.246      0.672  0.246   0
## ln_roads_km_c -0.063 0.358     -0.765   -0.063      0.638 -0.063   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                      mean       sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1284.692 14712.83      2.325   85.414   8040.646 3.540
## Phi for ui          0.343     0.28      0.009    0.269      0.927 0.012
## 
## Deviance Information Criterion (DIC) ...............: 62.65
## Deviance Information Criterion (DIC, saturated) ....: 45.93
## Effective number of parameters .....................: 3.17
## 
## Watanabe-Akaike information criterion (WAIC) ...: 62.90
## Effective number of parameters .................: 3.26
## 
## Marginal log-Likelihood:  1.57 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 14, Running = 0.403, Post = 0.124, Total = 14.5 
## Fixed effects:
##                                  mean     sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -56.307 30.490    -97.521  -65.611      6.254
## vandix_z                        0.069  0.370     -0.498    0.030      1.253
## ln_roads_km_c                  -0.993  1.583     -5.483   -0.811      2.554
## roads_prop_highway_arterial_z   0.654  0.868     -1.945    0.723      2.010
## ale_index_z                     2.985  3.395     -7.643    3.333      8.266
## canbics_index_z               -63.899 35.644   -113.806  -73.758      6.584
## population_100_c                0.758  0.628     -0.010    0.606      2.551
##                                  mode   kld
## (Intercept)                   -71.155 0.000
## vandix_z                       -0.003 0.005
## ln_roads_km_c                  -0.881 0.024
## roads_prop_highway_arterial_z   0.816 0.001
## ale_index_z                     3.650 0.003
## canbics_index_z               -88.299 0.000
## population_100_c                0.579 0.004
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                      mean       sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1381.138 16139.07      2.442   88.177   8570.257 3.794
## Phi for ui          0.344     0.28      0.009    0.269      0.927 0.012
## 
## Deviance Information Criterion (DIC) ...............: -1.29e+13
## Deviance Information Criterion (DIC, saturated) ....: -1.29e+13
## Effective number of parameters .....................: -1.29e+13
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1835.70
## Effective number of parameters .................: 889.55
## 
## Marginal log-Likelihood:  -7.91 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Model Set 3: Total Casualty Cyclist Claim Crashes

northwest_models_3 <- cma_models(northwest_sp@data,"Northwest","n_pedestrian_casualty_claims","vandix_z") 
## [1] "Fit unadjusted"
## [1] "Fit adjusted by road length"
## [1] "Fit adjusted by road length and covariates"
northwest_pedestrian <- clean_fixed_effects(northwest_models_3$fixed_effects)

map(northwest_models_3$models,~summary(.x))
## $Unadjusted
## Time used:
##     Pre = 14.3, Running = 0.387, Post = 0.118, Total = 14.8 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -1.469 0.344     -2.206   -1.450     -0.848 -1.400   0
## vandix_z     0.275 0.144     -0.006    0.274      0.563  0.274   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.544 0.898      0.516    1.324      3.903 0.993
## Phi for ui       0.341 0.234      0.027    0.295      0.849 0.080
## 
## Deviance Information Criterion (DIC) ...............: 142.59
## Deviance Information Criterion (DIC, saturated) ....: 90.73
## Effective number of parameters .....................: 19.88
## 
## Watanabe-Akaike information criterion (WAIC) ...: 147.89
## Effective number of parameters .................: 19.43
## 
## Marginal log-Likelihood:  -45.03 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted1
## Time used:
##     Pre = 13.9, Running = 0.382, Post = 0.137, Total = 14.4 
## Fixed effects:
##                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)   -1.681 0.402     -2.536   -1.661     -0.951 -1.602   0
## vandix_z       0.303 0.151      0.010    0.302      0.606  0.302   0
## ln_roads_km_c  0.252 0.242     -0.212    0.247      0.740  0.247   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                  mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 1.37 0.782      0.461    1.179      3.417 0.892
## Phi for ui       0.45 0.253      0.050    0.433      0.917 0.205
## 
## Deviance Information Criterion (DIC) ...............: 141.37
## Deviance Information Criterion (DIC, saturated) ....: 89.50
## Effective number of parameters .....................: 21.00
## 
## Watanabe-Akaike information criterion (WAIC) ...: 145.97
## Effective number of parameters .................: 19.60
## 
## Marginal log-Likelihood:  -49.40 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## 
## 
## $Adjusted2
## Time used:
##     Pre = 13.9, Running = 0.405, Post = 0.128, Total = 14.4 
## Fixed effects:
##                                 mean     sd 0.025quant 0.5quant 0.975quant
## (Intercept)                   -5.513  9.156    -23.468   -5.513     12.442
## vandix_z                       0.175  0.153     -0.126    0.174      0.478
## ln_roads_km_c                  0.143  0.441     -0.699    0.134      1.037
## roads_prop_highway_arterial_z  0.258  0.354     -0.452    0.263      0.941
## ale_index_z                    4.087  1.278      1.588    4.078      6.636
## canbics_index_z               -9.180 10.658    -30.075   -9.181     11.724
## population_100_c               0.349  0.162      0.032    0.347      0.674
##                                 mode kld
## (Intercept)                   -5.513   0
## vandix_z                       0.174   0
## ln_roads_km_c                  0.134   0
## roads_prop_highway_arterial_z  0.262   0
## ale_index_z                    4.078   0
## canbics_index_z               -9.181   0
## population_100_c               0.347   0
## 
## Random effects:
##   Name     Model
##     ui BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ui 2.332 1.734      0.617    1.858      6.956 1.247
## Phi for ui       0.305 0.218      0.024    0.255      0.804 0.071
## 
## Deviance Information Criterion (DIC) ...............: 159.78
## Deviance Information Criterion (DIC, saturated) ....: 107.92
## Effective number of parameters .....................: 28.84
## 
## Watanabe-Akaike information criterion (WAIC) ...: 189.93
## Effective number of parameters .................: 39.10
## 
## Marginal log-Likelihood:  -55.76 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
northwest_results <- bind_rows(northwest_all %>% mutate(Region = "Northwest",Outcome = "All Injury Claims"),
                         northwest_cyclist %>% mutate(Region = "Northwest",Outcome = "Cyclist Injury Claims"),
                         northwest_pedestrian %>% mutate(Region = "Northwest",Outcome = "Pedestrian Injury Claims")
                         ) %>%
  filter(variable == "vandix_z") %>%
  select(Region,Outcome,everything())

beepr::beep(8)

Results

Study Area

## Downloading: 12 kB     Downloading: 12 kB     Downloading: 12 kB     Downloading: 12 kB     Downloading: 16 kB     Downloading: 16 kB     Downloading: 16 kB     Downloading: 16 kB     Downloading: 24 kB     Downloading: 24 kB     Downloading: 24 kB     Downloading: 24 kB

Number of DAs

Total Injury Claims

Injury Claims Mean (SD)

Total Cycling Injury Claims

Cycling Injury Claims Mean (SD)

Total Pedestrian Injury Claims

Pedestrian Injury Claims Mean (SD)

Population

VANDIX Mean (SD)

Total Roads (km)

Highway/Arterial Roads Proportion Mean (SD)

No Highschool Prevalence Mean (SD)

Unemployment Rate Mean (SD)

Average Household Income Mean (SD)

Participation Rate Mean (SD)

University Degree Prevalence Mean (SD)

Lone Parent Family Prevalence Mean (SD)

Home Ownership Prevalence Mean (SD)

ALE Index Mean (SD)

CANBICS Index Mean (SD)

6,503

211,852

32.6 (70.8)

7,287

1.12 (2.77)

10,531

1.62 (3.46)

4,477,748

-0.016 (0.64)

47,865.21

0.146 (0.174)

0.158 (0.0805)

6.74 (4.62)

80100 (30400)

63.8 (10.9)

0.546 (0.118)

0.16 (0.0787)

0.691 (0.222)

0.692 (3.07)

3.15 (3.65)

The analysis included 6503 dissemination areas across 20 different metropolitan regions with a total of 211852 motor-vehicle crashes resulting in at least one injury, and 7287 such crashes involved at least one cyclist, while 10531 involved at least one pedestrian.

region_name

Number of DAs

Total Injury Claims

Injury Claims Mean (SD)

Total Cycling Injury Claims

Cycling Injury Claims Mean (SD)

Total Pedestrian Injury Claims

Pedestrian Injury Claims Mean (SD)

Population

VANDIX Mean (SD)

Total Roads (km)

Highway/Arterial Roads Proportion Mean (SD)

No Highschool Prevalence Mean (SD)

Unemployment Rate Mean (SD)

Average Household Income Mean (SD)

Participation Rate Mean (SD)

University Degree Prevalence Mean (SD)

Lone Parent Family Prevalence Mean (SD)

Home Ownership Prevalence Mean (SD)

ALE Index Mean (SD)

CANBICS Index Mean (SD)

Vancouver-Squamish

3,622

149,692

41.3 (84.8)

4,968

1.37 (3.24)

7,502

2.07 (4.02)

2,667,057

-0.14 (0.58)

15,056.825

0.159 (0.184)

0.143 (0.0767)

5.88 (3.58)

85400 (34200)

64.9 (10.1)

0.569 (0.121)

0.155 (0.0659)

0.671 (0.23)

1.63 (3.58)

4.57 (4.11)

Victoria

581

12,338

21.2 (34.2)

1,027

1.77 (3.17)

715

1.23 (2.54)

397,237

-0.24 (0.52)

3,261.164

0.116 (0.14)

0.121 (0.0581)

5.63 (3.87)

77300 (26900)

63.6 (11.1)

0.597 (0.0978)

0.153 (0.08)

0.654 (0.231)

0.607 (1.78)

2.12 (1.71)

Central Island-Powell River

611

10,762

17.6 (28.3)

271

0.444 (0.938)

547

0.895 (1.96)

357,193

0.2 (0.65)

5,652.122

0.13 (0.161)

0.177 (0.0756)

8.35 (5.15)

65700 (17200)

56.5 (11)

0.524 (0.0938)

0.167 (0.0948)

0.744 (0.188)

-1.07 (0.719)

0.588 (0.914)

Okanagan

466

13,236

28.4 (66.6)

485

1.04 (2.34)

564

1.21 (3.23)

336,628

0.12 (0.6)

4,759.053

0.13 (0.162)

0.164 (0.0651)

7.76 (4.73)

73000 (26800)

61.3 (12.9)

0.527 (0.0836)

0.158 (0.0831)

0.728 (0.208)

-0.843 (0.959)

1.9 (1.81)

Fraser Valley

472

15,730

33.3 (65.2)

360

0.763 (1.58)

725

1.54 (2.96)

309,493

0.31 (0.6)

3,678.908

0.0858 (0.138)

0.216 (0.0773)

7.05 (5.06)

75400 (21100)

63.6 (11.6)

0.444 (0.088)

0.166 (0.101)

0.741 (0.172)

-0.24 (1.73)

1.39 (1.3)

North Central

359

4,490

12.5 (25.2)

54

0.15 (0.544)

203

0.565 (1.62)

182,818

0.5 (0.71)

9,989.620

0.13 (0.145)

0.227 (0.0863)

10.5 (6.84)

77700 (21900)

67.9 (9.34)

0.453 (0.0895)

0.182 (0.1)

0.711 (0.236)

-1.32 (0.637)

0.582 (0.981)

Kamloops-Salmon Arm

204

4,069

19.9 (39.5)

89

0.436 (1.06)

183

0.897 (2.28)

133,847

0.18 (0.67)

2,744.521

0.245 (0.19)

0.17 (0.0717)

8.32 (5.27)

71500 (20400)

62.3 (9.95)

0.511 (0.0944)

0.178 (0.101)

0.74 (0.201)

-0.714 (1.22)

1.4 (1.47)

Southeast

110

1,037

9.43 (17.8)

24

0.218 (0.612)

51

0.464 (1.23)

60,427

0.079 (0.58)

1,598.577

0.128 (0.156)

0.156 (0.0625)

8.06 (4.35)

68700 (16800)

62.2 (8.54)

0.549 (0.0921)

0.165 (0.0835)

0.743 (0.203)

-1.24 (0.743)

0.321 (0.535)

Northwest

78

498

6.38 (11.8)

9

0.115 (0.36)

41

0.526 (1.05)

33,048

0.68 (0.91)

1,124.422

0.156 (0.223)

0.247 (0.098)

12.5 (9.31)

74400 (17000)

67.8 (7.67)

0.457 (0.102)

0.196 (0.0949)

0.7 (0.189)

-1.49 (0.483)

0 (0)

Association between Deprivation Index and Traffic Injury

Crude Rates

All Motor Vehicle Traffic Injury Crashes

Bicycle-Motor Vehicle Traffic Injury Crashes

Pedestrian-Motor Vehicle Traffic Injury Crashes

Spatial Models

All Motor Vehicle Traffic Injury Crashes

Relationship between VanDIX score and total number of motor-vehicle crash insurance claims from 2019-2023.

Variable

Region

Unadjusted IRR (95% CI)

Fully Adjusted IRR (95% CI)a

VanDIX (SD)

Vancouver-Squamish

1.24 (1.16, 1.32)

1.25 (1.2, 1.3)

Victoria

1.22 (1.07, 1.38)

1.22 (1.11, 1.35)

Central Island-Powell River

1.34 (1.19, 1.52)

1.37 (1.26, 1.5)

Okanagan

1.41 (1.18, 1.69)

1.43 (1.26, 1.62)

Fraser Valley

1.26 (1.07, 1.49)

1.15 (1.02, 1.29)

North Central

1.21 (1.04, 1.41)

1.22 (1.07, 1.39)

Kamloops-Salmon Arm

1.28 (1.02, 1.61)

1.23 (1.04, 1.46)

Southeast

1.42 (1.06, 1.91)

1.42 (1.1, 1.84)

Northwest

1.32 (1.04, 1.68)

1.25 (1.01, 1.54)

aAdjusted for (i) total length of roads; (ii) proportion of roads classified as highway or arterial; (iii) the Canadian Active Living Environment Index; (iv) the Canadian Bikeway Comfort and Safety Index; (v) Total population

Bicycle-Motor Vehicle Traffic Injury Crashes

Relationship between VanDIX score and total number of motor-vehicle - cyclist crash insurance claims from 2019-2023.

Variable

Region

Unadjusted IRR (95% CI)

Fully Adjusted IRR (95% CI)a

VanDIX (SD)

Vancouver-Squamish

1.16 (1.08, 1.25)

1.2 (1.13, 1.28)

Victoria

1.16 (0.983, 1.37)

1.13 (0.962, 1.32)

Central Island-Powell River

1.4 (1.21, 1.62)

1.42 (1.22, 1.66)

Okanagan

1.9 (1.54, 2.36)

1.76 (1.45, 2.15)

Fraser Valley

1.4 (1.12, 1.75)

1.08 (0.87, 1.35)

North Central

1.36 (1.04, 1.79)

1.43 (1.07, 1.92)

Kamloops-Salmon Arm

1.45 (1.07, 1.95)

1.48 (1.03, 2.07)

Southeast

1.82 (1.17, 2.89)

2.29 (1.4, 3.76)

Northwest

1.27 (0.841, 1.91)

1.07 (0.607, 3.5)

aAdjusted for (i) total length of roads; (ii) proportion of roads classified as highway or arterial; (iii) the Canadian Active Living Environment Index; (iv) the Canadian Bikeway Comfort and Safety Index; (v) Total population

Pedestrian-Motor Vehicle Traffic Injury Crashes

Table 3: Relationship between VanDIX score and total number of motor-vehicle - pedestrian crash insurance claims from 2019-2023.

Variable

Region

Unadjusted IRR (95% CI)

Fully Adjusted IRR (95% CI)a

VanDIX (SD)

Vancouver-Squamish

1.3 (1.21, 1.4)

1.31 (1.23, 1.39)

Victoria

1.42 (1.2, 1.67)

1.37 (1.17, 1.61)

Central Island-Powell River

1.74 (1.5, 2.03)

1.68 (1.45, 1.95)

Okanagan

1.95 (1.57, 2.42)

1.75 (1.4, 2.18)

Fraser Valley

1.65 (1.37, 1.99)

1.32 (1.11, 1.57)

North Central

1.47 (1.19, 1.8)

1.4 (1.13, 1.73)

Kamloops-Salmon Arm

1.92 (1.41, 2.64)

1.94 (1.44, 2.63)

Southeast

1.72 (1.08, 2.78)

1.83 (1.12, 3.02)

Northwest

1.32 (0.994, 1.76)

1.19 (0.882, 1.61)

aAdjusted for (i) total length of roads; (ii) proportion of roads classified as highway or arterial; (iii) the Canadian Active Living Environment Index; (iv) the Canadian Bikeway Comfort and Safety Index; (v) Total population

Supplementary Material

Built Environment by VanDIX

  • Of the 9 regions, we observed 5 had on average a positive relationship with deprivation index and the proportion of roads in the neighborhood that were arterial or highway

  • Of the 9 regions, we observed 8 had on average a positive relationship with deprivation index and thactive living environment index score

  • Of the 9 regions, we observed 8 had on average a positive relationship with deprivation index and increased access to better bicycle facilities.

Model Details

Vancouver Model Details

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

2.800 (2.770, 2.830)

3.250 (3.220, 3.280)

2.990 (2.960, 3.020)

vandix_z

0.211 (0.147, 0.276)

0.327 (0.276, 0.378)

0.222 (0.180, 0.265)

ln_roads_km_c

1.290 (1.240, 1.350)

0.990 (0.937, 1.040)

roads_prop_highway_arterial_z

0.593 (0.562, 0.624)

ale_index_z

0.109 (0.036, 0.182)

canbics_index_z

0.209 (0.137, 0.282)

population_100_c

0.017 (0.012, 0.022)

Hyper Parameters

Precision for ui

0.311 (0.284, 0.339)

0.508 (0.470, 0.548)

0.838 (0.765, 0.915)

Phi for ui

0.794 (0.751, 0.833)

0.811 (0.778, 0.841)

0.777 (0.733, 0.818)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-0.523 (-0.588, -0.461)

-0.134 (-0.186, -0.083)

-0.291 (-0.352, -0.231)

vandix_z

0.148 (0.075, 0.221)

0.217 (0.154, 0.280)

0.185 (0.123, 0.248)

ln_roads_km_c

1.070 (1.000, 1.140)

0.844 (0.767, 0.921)

roads_prop_highway_arterial_z

0.301 (0.254, 0.348)

ale_index_z

0.041 (-0.048, 0.130)

canbics_index_z

0.092 (-0.001, 0.186)

population_100_c

0.022 (0.016, 0.028)

Hyper Parameters

Precision for ui

0.567 (0.489, 0.653)

0.622 (0.532, 0.718)

0.822 (0.692, 0.964)

Phi for ui

0.633 (0.534, 0.724)

0.934 (0.874, 0.974)

0.871 (0.790, 0.931)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-0.035 (-0.089, 0.018)

0.319 (0.275, 0.363)

0.115 (0.060, 0.168)

vandix_z

0.262 (0.191, 0.333)

0.304 (0.242, 0.367)

0.268 (0.208, 0.328)

ln_roads_km_c

0.997 (0.930, 1.060)

0.717 (0.641, 0.793)

roads_prop_highway_arterial_z

0.394 (0.351, 0.438)

ale_index_z

0.232 (0.141, 0.322)

canbics_index_z

0.042 (-0.051, 0.135)

population_100_c

0.027 (0.021, 0.034)

Hyper Parameters

Precision for ui

0.496 (0.430, 0.567)

0.450 (0.396, 0.506)

0.713 (0.602, 0.835)

Phi for ui

0.663 (0.573, 0.746)

0.950 (0.907, 0.979)

0.850 (0.763, 0.917)

Victoria Model Details

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

2.340 (2.270, 2.410)

2.360 (2.310, 2.410)

2.510 (2.420, 2.600)

vandix_z

0.196 (0.068, 0.325)

0.326 (0.222, 0.430)

0.202 (0.104, 0.301)

ln_roads_km_c

1.240 (1.100, 1.380)

1.010 (0.850, 1.160)

roads_prop_highway_arterial_z

0.462 (0.367, 0.558)

ale_index_z

0.409 (0.081, 0.739)

canbics_index_z

0.340 (0.032, 0.646)

population_100_c

0.015 (-0.002, 0.031)

Hyper Parameters

Precision for ui

0.490 (0.390, 0.599)

0.642 (0.547, 0.744)

0.845 (0.679, 1.020)

Phi for ui

0.857 (0.721, 0.947)

0.981 (0.930, 0.999)

0.948 (0.838, 0.994)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-0.141 (-0.282, -0.006)

-0.144 (-0.273, -0.018)

0.039 (-0.120, 0.193)

vandix_z

0.150 (-0.017, 0.318)

0.250 (0.094, 0.406)

0.118 (-0.038, 0.275)

ln_roads_km_c

1.100 (0.871, 1.320)

0.887 (0.635, 1.140)

roads_prop_highway_arterial_z

0.308 (0.165, 0.451)

ale_index_z

0.401 (-0.033, 0.838)

canbics_index_z

0.587 (0.169, 1.010)

population_100_c

0.030 (0.004, 0.058)

Hyper Parameters

Precision for ui

0.679 (0.511, 0.882)

0.663 (0.489, 0.868)

0.867 (0.622, 1.160)

Phi for ui

0.623 (0.414, 0.806)

0.902 (0.739, 0.984)

0.828 (0.616, 0.957)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-0.439 (-0.589, -0.296)

-0.430 (-0.573, -0.292)

-0.313 (-0.493, -0.141)

vandix_z

0.350 (0.184, 0.516)

0.452 (0.292, 0.612)

0.317 (0.156, 0.478)

ln_roads_km_c

0.854 (0.625, 1.090)

0.578 (0.313, 0.845)

roads_prop_highway_arterial_z

0.384 (0.228, 0.539)

ale_index_z

0.540 (0.080, 1.010)

canbics_index_z

0.326 (-0.136, 0.786)

population_100_c

0.037 (0.010, 0.064)

Hyper Parameters

Precision for ui

0.633 (0.458, 0.847)

0.619 (0.458, 0.816)

0.758 (0.543, 1.020)

Phi for ui

0.791 (0.569, 0.937)

0.950 (0.822, 0.996)

0.924 (0.745, 0.993)

Central Island - Powell River Model Details

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

1.930 (1.830, 2.020)

1.360 (1.240, 1.470)

1.960 (1.730, 2.190)

vandix_z

0.295 (0.175, 0.416)

0.435 (0.330, 0.539)

0.317 (0.227, 0.408)

ln_roads_km_c

1.020 (0.888, 1.160)

0.650 (0.506, 0.794)

roads_prop_highway_arterial_z

0.668 (0.575, 0.761)

ale_index_z

1.420 (0.937, 1.900)

canbics_index_z

-0.780 (-1.230, -0.328)

population_100_c

0.063 (0.037, 0.090)

Hyper Parameters

Precision for ui

0.410 (0.309, 0.537)

0.442 (0.345, 0.558)

0.670 (0.517, 0.852)

Phi for ui

0.571 (0.383, 0.741)

0.736 (0.613, 0.836)

0.726 (0.598, 0.832)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.430 (-1.650, -1.220)

-1.730 (-2.010, -1.460)

-1.150 (-1.570, -0.753)

vandix_z

0.335 (0.188, 0.483)

0.399 (0.250, 0.548)

0.354 (0.203, 0.504)

ln_roads_km_c

0.506 (0.280, 0.737)

0.292 (0.018, 0.568)

roads_prop_highway_arterial_z

0.431 (0.272, 0.592)

ale_index_z

1.160 (0.350, 1.970)

canbics_index_z

-0.321 (-1.040, 0.394)

population_100_c

0.060 (0.014, 0.107)

Hyper Parameters

Precision for ui

1.060 (0.680, 1.600)

0.927 (0.561, 1.460)

1.180 (0.691, 1.910)

Phi for ui

0.320 (0.089, 0.626)

0.551 (0.249, 0.823)

0.496 (0.197, 0.788)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.250 (-1.490, -1.040)

-1.600 (-1.880, -1.330)

-0.755 (-1.170, -0.361)

vandix_z

0.556 (0.405, 0.710)

0.611 (0.456, 0.769)

0.518 (0.370, 0.666)

ln_roads_km_c

0.653 (0.401, 0.908)

0.328 (0.052, 0.605)

roads_prop_highway_arterial_z

0.552 (0.394, 0.713)

ale_index_z

2.370 (1.610, 3.150)

canbics_index_z

-1.150 (-1.840, -0.468)

population_100_c

0.115 (0.071, 0.160)

Hyper Parameters

Precision for ui

0.605 (0.439, 0.830)

0.468 (0.302, 0.701)

0.749 (0.489, 1.120)

Phi for ui

0.206 (0.030, 0.511)

0.587 (0.295, 0.830)

0.446 (0.166, 0.736)

Okanagan Model Details

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

2.160 (2.070, 2.260)

1.580 (1.480, 1.680)

1.960 (1.800, 2.130)

vandix_z

0.344 (0.163, 0.525)

0.464 (0.318, 0.610)

0.355 (0.228, 0.484)

ln_roads_km_c

1.150 (1.000, 1.290)

0.949 (0.801, 1.100)

roads_prop_highway_arterial_z

0.605 (0.499, 0.711)

ale_index_z

0.510 (-0.006, 1.030)

canbics_index_z

0.016 (-0.406, 0.438)

population_100_c

0.031 (0.014, 0.048)

Hyper Parameters

Precision for ui

0.237 (0.172, 0.317)

0.241 (0.200, 0.286)

0.339 (0.271, 0.412)

Phi for ui

0.843 (0.704, 0.934)

0.985 (0.950, 0.998)

0.979 (0.929, 0.998)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.110 (-1.340, -0.903)

-1.510 (-1.790, -1.240)

-1.010 (-1.350, -0.697)

vandix_z

0.642 (0.429, 0.858)

0.705 (0.500, 0.913)

0.566 (0.369, 0.765)

ln_roads_km_c

0.721 (0.484, 0.961)

0.616 (0.348, 0.887)

roads_prop_highway_arterial_z

0.395 (0.221, 0.569)

ale_index_z

0.703 (-0.040, 1.450)

canbics_index_z

0.354 (-0.220, 0.925)

population_100_c

0.036 (0.008, 0.065)

Hyper Parameters

Precision for ui

0.426 (0.269, 0.634)

0.339 (0.230, 0.474)

0.428 (0.288, 0.608)

Phi for ui

0.788 (0.568, 0.933)

0.951 (0.840, 0.995)

0.961 (0.858, 0.997)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.040 (-1.280, -0.822)

-1.440 (-1.720, -1.170)

-0.861 (-1.180, -0.554)

vandix_z

0.668 (0.452, 0.884)

0.731 (0.497, 0.965)

0.558 (0.338, 0.780)

ln_roads_km_c

0.805 (0.541, 1.070)

0.632 (0.337, 0.927)

roads_prop_highway_arterial_z

0.575 (0.404, 0.748)

ale_index_z

1.290 (0.531, 2.070)

canbics_index_z

-0.245 (-0.877, 0.384)

population_100_c

0.049 (0.019, 0.081)

Hyper Parameters

Precision for ui

0.525 (0.374, 0.720)

0.372 (0.232, 0.561)

0.508 (0.310, 0.786)

Phi for ui

0.328 (0.121, 0.593)

0.751 (0.489, 0.925)

0.707 (0.408, 0.911)

Fraser Valley Model Details

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

2.480 (2.400, 2.560)

2.270 (2.180, 2.350)

2.720 (2.600, 2.840)

vandix_z

0.232 (0.065, 0.399)

0.285 (0.156, 0.415)

0.140 (0.022, 0.257)

ln_roads_km_c

1.100 (0.972, 1.240)

0.862 (0.721, 1.000)

roads_prop_highway_arterial_z

0.580 (0.464, 0.695)

ale_index_z

0.245 (-0.004, 0.494)

canbics_index_z

0.183 (-0.097, 0.458)

population_100_c

0.027 (0.003, 0.051)

Hyper Parameters

Precision for ui

0.299 (0.258, 0.344)

0.712 (0.585, 0.856)

0.953 (0.760, 1.170)

Phi for ui

0.767 (0.001, 1.000)

0.802 (0.680, 0.893)

0.805 (0.658, 0.913)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.220 (-1.480, -0.984)

-1.360 (-1.610, -1.120)

-0.791 (-1.080, -0.513)

vandix_z

0.338 (0.113, 0.557)

0.286 (0.063, 0.509)

0.081 (-0.139, 0.299)

ln_roads_km_c

0.753 (0.505, 1.000)

0.467 (0.207, 0.731)

roads_prop_highway_arterial_z

0.601 (0.404, 0.800)

ale_index_z

0.695 (0.276, 1.110)

canbics_index_z

0.091 (-0.494, 0.678)

population_100_c

0.046 (0.001, 0.091)

Hyper Parameters

Precision for ui

0.679 (0.485, 0.930)

0.623 (0.427, 0.874)

0.897 (0.588, 1.320)

Phi for ui

0.315 (0.098, 0.593)

0.752 (0.479, 0.934)

0.611 (0.278, 0.882)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-0.585 (-0.779, -0.402)

-0.730 (-0.915, -0.552)

-0.210 (-0.437, 0.010)

vandix_z

0.503 (0.318, 0.689)

0.447 (0.269, 0.627)

0.279 (0.109, 0.450)

ln_roads_km_c

0.813 (0.608, 1.020)

0.421 (0.199, 0.646)

roads_prop_highway_arterial_z

0.619 (0.459, 0.781)

ale_index_z

0.488 (0.140, 0.837)

canbics_index_z

0.025 (-0.445, 0.494)

population_100_c

0.059 (0.023, 0.097)

Hyper Parameters

Precision for ui

0.740 (0.558, 0.965)

0.660 (0.476, 0.882)

1.040 (0.732, 1.440)

Phi for ui

0.372 (0.166, 0.610)

0.852 (0.638, 0.970)

0.675 (0.392, 0.892)

North Central Model Details

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

1.320 (1.160, 1.480)

0.531 (0.261, 0.795)

0.969 (0.535, 1.400)

vandix_z

0.191 (0.042, 0.340)

0.272 (0.132, 0.413)

0.201 (0.071, 0.332)

ln_roads_km_c

0.698 (0.521, 0.876)

0.542 (0.371, 0.713)

roads_prop_highway_arterial_z

0.616 (0.447, 0.784)

ale_index_z

-0.263 (-0.837, 0.313)

canbics_index_z

0.384 (-0.228, 0.994)

population_100_c

0.034 (-0.008, 0.076)

Hyper Parameters

Precision for ui

0.296 (0.218, 0.391)

0.229 (0.172, 0.294)

0.269 (0.204, 0.341)

Phi for ui

0.788 (0.650, 0.890)

0.950 (0.871, 0.990)

0.962 (0.890, 0.994)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-3.130 (-3.910, -2.500)

-3.340 (-4.190, -2.590)

-2.520 (-3.540, -1.580)

vandix_z

0.311 (0.043, 0.584)

0.350 (0.069, 0.637)

0.356 (0.068, 0.651)

ln_roads_km_c

0.160 (-0.154, 0.483)

0.211 (-0.141, 0.570)

roads_prop_highway_arterial_z

0.504 (0.115, 0.898)

ale_index_z

-0.474 (-1.870, 0.929)

canbics_index_z

2.030 (0.911, 3.170)

population_100_c

-0.047 (-0.197, 0.104)

Hyper Parameters

Precision for ui

0.649 (0.347, 1.130)

0.637 (0.340, 1.110)

0.833 (0.399, 1.590)

Phi for ui

0.158 (0.024, 0.438)

0.199 (0.032, 0.522)

0.306 (0.063, 0.682)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.930 (-2.340, -1.570)

-2.060 (-2.570, -1.600)

-2.110 (-3.010, -1.270)

vandix_z

0.382 (0.177, 0.589)

0.392 (0.182, 0.604)

0.335 (0.126, 0.548)

ln_roads_km_c

0.120 (-0.147, 0.409)

-0.063 (-0.372, 0.265)

roads_prop_highway_arterial_z

0.650 (0.364, 0.939)

ale_index_z

-1.240 (-2.500, -0.021)

canbics_index_z

0.642 (-0.340, 1.640)

population_100_c

0.022 (-0.064, 0.109)

Hyper Parameters

Precision for ui

0.548 (0.346, 0.833)

0.512 (0.302, 0.818)

0.522 (0.291, 0.864)

Phi for ui

0.359 (0.106, 0.681)

0.451 (0.131, 0.801)

0.560 (0.204, 0.875)

Kamloops-Salmon Arm

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

1.950 (1.790, 2.110)

1.430 (1.230, 1.620)

1.130 (0.898, 1.370)

vandix_z

0.250 (0.024, 0.476)

0.258 (0.055, 0.461)

0.208 (0.042, 0.375)

ln_roads_km_c

0.850 (0.644, 1.060)

0.491 (0.302, 0.681)

roads_prop_highway_arterial_z

0.534 (0.424, 0.645)

ale_index_z

0.091 (-0.413, 0.594)

canbics_index_z

-0.712 (-1.200, -0.224)

population_100_c

0.099 (0.062, 0.136)

Hyper Parameters

Precision for ui

0.647 (0.498, 0.825)

0.776 (0.590, 0.999)

1.180 (0.871, 1.540)

Phi for ui

0.421 (0.223, 0.639)

0.625 (0.425, 0.802)

0.819 (0.620, 0.948)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.620 (-2.050, -1.230)

-2.040 (-2.620, -1.520)

-1.920 (-2.590, -1.320)

vandix_z

0.372 (0.064, 0.668)

0.394 (0.069, 0.709)

0.391 (0.034, 0.726)

ln_roads_km_c

0.586 (0.200, 0.990)

0.195 (-0.224, 0.624)

roads_prop_highway_arterial_z

0.395 (0.125, 0.668)

ale_index_z

0.630 (-0.340, 1.580)

canbics_index_z

-0.284 (-1.240, 0.683)

population_100_c

0.141 (0.067, 0.217)

Hyper Parameters

Precision for ui

1.060 (0.563, 1.860)

1.070 (0.551, 1.910)

1.850 (0.722, 4.150)

Phi for ui

0.209 (0.034, 0.528)

0.519 (0.160, 0.871)

0.619 (0.170, 0.950)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.460 (-1.890, -1.090)

-1.850 (-2.350, -1.390)

-1.580 (-2.170, -1.040)

vandix_z

0.654 (0.343, 0.970)

0.693 (0.377, 1.010)

0.665 (0.362, 0.966)

ln_roads_km_c

0.620 (0.269, 0.980)

0.248 (-0.120, 0.622)

roads_prop_highway_arterial_z

0.389 (0.148, 0.633)

ale_index_z

1.570 (0.675, 2.450)

canbics_index_z

-0.983 (-1.880, -0.103)

population_100_c

0.111 (0.039, 0.184)

Hyper Parameters

Precision for ui

0.626 (0.404, 0.925)

0.676 (0.426, 1.020)

1.080 (0.602, 1.820)

Phi for ui

0.267 (0.075, 0.549)

0.474 (0.190, 0.778)

0.373 (0.064, 0.783)

Southeast

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

1.270 (1.050, 1.480)

0.651 (0.321, 0.973)

1.560 (0.776, 2.340)

vandix_z

0.352 (0.057, 0.648)

0.386 (0.111, 0.664)

0.350 (0.094, 0.609)

ln_roads_km_c

0.779 (0.474, 1.090)

0.490 (0.142, 0.844)

roads_prop_highway_arterial_z

0.627 (0.331, 0.921)

ale_index_z

1.630 (0.032, 3.250)

canbics_index_z

-0.712 (-2.390, 0.955)

population_100_c

0.152 (0.040, 0.265)

Hyper Parameters

Precision for ui

0.750 (0.528, 1.030)

0.827 (0.571, 1.150)

0.994 (0.677, 1.400)

Phi for ui

0.396 (0.151, 0.689)

0.691 (0.377, 0.923)

0.771 (0.457, 0.961)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-2.000 (-2.780, -1.390)

-2.450 (-3.360, -1.640)

-0.882 (-2.090, 0.301)

vandix_z

0.597 (0.156, 1.060)

0.755 (0.271, 1.250)

0.828 (0.337, 1.330)

ln_roads_km_c

0.468 (-0.008, 0.948)

0.202 (-0.535, 0.939)

roads_prop_highway_arterial_z

0.331 (-0.323, 0.989)

ale_index_z

1.940 (-0.715, 4.620)

canbics_index_z

0.074 (-2.440, 2.590)

population_100_c

0.288 (0.105, 0.473)

Hyper Parameters

Precision for ui

11.100 (0.759, 67.100)

60.400 (0.709, 408.000)

5,990.000 (0.933, 22,800.000)

Phi for ui

0.277 (0.018, 0.796)

0.357 (0.028, 0.880)

0.353 (0.018, 0.906)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.870 (-2.560, -1.290)

-2.100 (-2.900, -1.400)

-1.380 (-3.000, 0.125)

vandix_z

0.545 (0.075, 1.020)

0.597 (0.113, 1.080)

0.606 (0.111, 1.110)

ln_roads_km_c

0.295 (-0.211, 0.801)

-0.070 (-0.840, 0.685)

roads_prop_highway_arterial_z

0.307 (-0.349, 0.958)

ale_index_z

0.615 (-2.210, 3.390)

canbics_index_z

-0.167 (-2.980, 2.700)

population_100_c

0.233 (0.010, 0.464)

Hyper Parameters

Precision for ui

0.728 (0.368, 1.310)

0.753 (0.372, 1.380)

0.736 (0.349, 1.390)

Phi for ui

0.165 (0.011, 0.546)

0.204 (0.014, 0.634)

0.227 (0.016, 0.679)

Northwest

All

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

0.682 (0.281, 1.050)

0.162 (-0.275, 0.573)

-5.270 (-21.900, 11.400)

vandix_z

0.281 (0.042, 0.521)

0.362 (0.138, 0.591)

0.222 (0.013, 0.433)

ln_roads_km_c

0.780 (0.469, 1.100)

0.423 (-0.130, 0.986)

roads_prop_highway_arterial_z

0.228 (-0.209, 0.663)

ale_index_z

3.510 (1.690, 5.290)

canbics_index_z

-10.800 (-30.200, 8.550)

population_100_c

0.372 (0.130, 0.614)

Hyper Parameters

Precision for ui

0.632 (0.381, 0.999)

0.559 (0.323, 0.885)

1.160 (0.688, 1.850)

Phi for ui

0.367 (0.054, 0.785)

0.810 (0.430, 0.983)

0.192 (0.002, 0.776)

Cyclists

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-2.580 (-3.500, -1.670)

-2.630 (-3.650, -1.610)

-56.300 (-97.500, 6.250)

vandix_z

0.237 (-0.173, 0.649)

0.246 (-0.180, 0.672)

0.069 (-0.498, 1.250)

ln_roads_km_c

-0.063 (-0.765, 0.638)

-0.993 (-5.480, 2.550)

roads_prop_highway_arterial_z

0.654 (-1.950, 2.010)

ale_index_z

2.980 (-7.640, 8.270)

canbics_index_z

-63.900 (-114.000, 6.580)

population_100_c

0.758 (-0.010, 2.550)

Hyper Parameters

Precision for ui

1,270.000 (2.360, 7,980.000)

1,280.000 (2.330, 8,040.000)

1,380.000 (2.440, 8,570.000)

Phi for ui

0.343 (0.009, 0.927)

0.343 (0.009, 0.927)

0.344 (0.009, 0.927)

Pedestrians

term

Unadjusted

Adjusted1

Adjusted2

Fixed Effects

(Intercept)

-1.470 (-2.210, -0.848)

-1.680 (-2.540, -0.951)

-5.510 (-23.500, 12.400)

vandix_z

0.275 (-0.006, 0.563)

0.303 (0.010, 0.606)

0.175 (-0.126, 0.478)

ln_roads_km_c

0.252 (-0.212, 0.740)

0.143 (-0.699, 1.040)

roads_prop_highway_arterial_z

0.258 (-0.452, 0.941)

ale_index_z

4.090 (1.590, 6.640)

canbics_index_z

-9.180 (-30.100, 11.700)

population_100_c

0.349 (0.032, 0.674)

Hyper Parameters

Precision for ui

1.540 (0.516, 3.900)

1.370 (0.461, 3.420)

2.330 (0.617, 6.960)

Phi for ui

0.341 (0.027, 0.849)

0.450 (0.050, 0.917)

0.305 (0.024, 0.804)

end_time <- Sys.time()
end_time - start_time
## Time difference of 23.88309 mins